Playing with the Chaos Game: non-Euclidean Sierpinski attractor in spherical coordinates

In my former post about The Chaos Game I did a non-Euclidean Sierpinski attractor by using polar coordinates in the XY plane (2D). Now I have tried to make a new "flavor" by using spherical coordinates in the XYZ region (3D). 

In this case, I have divided the sphere into 8 quarters, and the Sierpinski attractor is generated in the surface of a quarter of the sphere and the rules of the Chaos Game regarding how to manage the angles are exactly the same rules as for the polar coordinates version of my former post, but with one more dimension (in other words, one more angle to play with). Specifically the three attractor points are at A=(phi=0,theta=0,r) (east of the XY plane), B=(phi=3PI/2,theta=0,r) (south of the XY plane) and C=(phi=7PI/4,tetha=PI,r) (north of the XZ plane). This produces a non-Euclidean spherical Sierpinski attractor in this fashion:

 

The next step is just replicating the pattern symmetrically on the surfaces of the rest of quarters of the sphere, and this is the result! (it is an animated gif, it could take a little bit to load):

Two more views:

I am really very happy with the results! The spherical symmetry is beautiful. It is easy to do bad calculations when working with polar and spherical coordinates. If possible I will try to add more "flavors" of The Chaos Game. If somebody is interested in the Python code, please just let me know!

Hobbymath's wanderings #9

A random walk through fractal dimensions.

Wikipedia: the Chaos Game.

Wolfram: the Chaos Game.

About the Chaos Game.

Examples of Iterated Function Systems.

The Beauty of Roots (polynomials of degree ≤ 5).

Barnsley Fern.

Superfractals: the Chaos Game.

Adding a little more chaos to the Chaos Game.

The Chaos Game in polar coordinates: a nice non-Euclidean Sierpinski attractor

I am learning the basic concepts about the Chaos Game (I did a previous question about the same topic at MSE here), the method to create fractals elaborated by professor Michael Barnsley.

The basic example (euclidean-cartesian) requires three attractor vertices, ABC, and as Wikipedia explains: "The fractal is created by iteratively creating a sequence of points, starting with the initial random point, in which each point in the sequence is a given fraction of the distance between the previous point and one of the vertices of the polygon; the vertex is chosen at random in each iteration."

The sequence of points (x,y) will produce the pattern of the Sierpinski attractor:

When I saw that construction, I wondered if there was another way of obtaining an attractor with two variables and three attractor points, and thought about polar coordinates. I did not find any references, so I did my own version of the Chaos Game for polar coordinates.

In this version, the points are (theta,r) (the angle in radians and the radius). And the three attractor points are A=(0,0), B=(5PI/4,1) and C=(7PI/4,10^4).

This is the algorithm I have implemented:

 And this is the result for 10^7 iterations, a nice Sierpinski attractor:

The Sierpinski attractor appears, but there is a distortion due to the use of polar coordinates. In the case of the polar coordinates the point (theta, r) has two options to arrive to the attractor point: using the clockwise angle or the counterclockwise angle. Only selecting the shortest path (in this case the smallest angle of both options) the Sierpinski attractor is shown. It looked like a (sorry if I am not using the appropriated words) "non-Euclidean" Sierpinski attractor.

This is another example locating the attractor points in the same axis: A=(0,0), B=(PI-PI/8,1) and C=(7PI/4,10^4).

I wondered if, according to the results, can be said that the Sierpinski attractor in the polar coordinates pattern is a "non-Euclidean" version of the attractor (or it is only a "Euclidean" distortion of the original attractor). Fortunately a fellow from MSE confirmed that it is fair to say that is is a non-Euclidean version! Here is the link to my question and the answer. I am very happy because I think this might be the first version of the Chaos Game in polar coordinates showing the Sierpinski attractor in a non-Euclidean "flavor", I could be wrong but initially I did not find it in Internet!

An interesting sequence obtained from an orbits simulator: Fibonacci (II)

In the previous post I wrote an orbit simulator that produces a very intriguing sequence. I am still reviewing it, and still did not see any special property of it, but using the same orbits algorithm and instead of  using the sequence of natural numbers, using Fibonacci numbers the results are much more interesting!

Here are the results of the first terms, when the sequence is the Fibonacci numbers, due to the relationship of the Fibonacci terms, only two orbits are generated. A third orbit will never be possible!

Observing the two orbits it is easy to see that the odd terms of the first orbit divide the odd terms of the second orbit in the same position. And here is the surprise, when dividing the second orbit terms by the first orbit terms, the appearing sequence is:

{1, 7, 29, 123, 521, 2207, 9349, 39603, 167761, 710647, 3010349, 12752043, 54018521, ...}

Which is the generalized Pell equation with second term of 7, and can be found at OEIS here.

In the other hand the odd terms of the second orbit are Fibonacci(6n+5).

The second orbit appears as well as a sequence at OEIS (A015448). a(0)=1, a(1)=1, and a(n) = 4*a(n-1) + a(n-2) for n>=2.

About the first orbit, only the odd terms appear as A033887:  a(n) = Fibonacci(3n+1).

The results are quite interesting and I hope they will help me to find some relationship when using the natural numbers sequence instead of Fibonacci.

An interesting sequence obtained from an orbits simulator

I have obtained an interesting sequence by making a simulation of a system with a total of k "orbits", system_k = {o_0,o_1..o_k-1} where the orbits are filled with natural numbers following the algorithm explained below. The orbit closer to the "nucleus" is o_0. Applying the rules below I can obtain a sequence that I suspected would be easy to find at OEIS, but it was not there.

The main rule to generate the sequence is: an orbit o_i can never have a set of numbers whose sum is greater than o_i-1 (in other words, the more internal orbit it is, the "heavier" it is) and the most external orbit only can have one number. The rules to settle numbers follow this algorithm:

1. number to settle n = 1.

2. current orbit = the most external orbit = o_k-1

3. If the current orbit is not o_0 and the sum of the numbers currently settled in the current orbit plus n are less or equal than the sum of the numbers currently settled in the immediately lower closer orbit then add n to the current orbit. If the current orbit is the most external one, finish.

4. In other case, if the current orbit is not o_0 then current orbit = immediately lower closer orbit and try again step (3). If the current orbit is o_0 directly settle n into it. If it was a system of total orbits = 1 then finish. If not, n=n+1 and go to step (2).

By doing this the algorithm finishes when a number n is settled for the first time the most external orbit of the system.

(Q) What I wanted to simulate is: for the interval of [1,T] orbits, how many numbers are necessary to settle until one number is finally settled in the outer orbit? 

E.g. this are the results applying the algorithm to some orbits:

It is easy to see that the system_i is a subset of the system_i+1. It means that finding the answer to the question (Q) is the same as generating a T orbits system (system_t) and checking which one was the first number n settled on each orbit of the system.

When the test is run this is the sequence S of numbers required to settle a number in the outer orbit for a number of [1..49] orbits:

S={1,3,7,11,16,25,32,40,49,59,73,85,98,112,127,147,164,182,201,221,242,272,295,319,344,370,397,425,454,484,515,553,586,620,655,691}

The sequence can not be found at OEIS, but the idea is simple so I suspect that probably there is this same sequence somewhere and can not find it.

I wonder if there is a reference about a similar sequence like this, applying a similar algorithm, and if the sequence S is somehow predictable (there is a closed-form expression).

Some extra information:

This is how system_31 loos like:

And this is the Python code in case somebody wanted to play with it:

def orbits():
    from sympy import totient,mobius
    from gmpy2 import is_prime,is_square
  
    orbits=[[],[]]
    total_orbits = 2
    current_orbit = total_orbits
    for n in range (1,100):
        while True:
            valn = n
            if (current_orbit>1) and (sum(orbits[current_orbit-1])+valn) <= (sum(orbits[current_orbit-2])):
                orbits[current_orbit-1].append(valn)
                break

            else:
                current_orbit = current_orbit -1
                if current_orbit == 1 or current_orbit==0:
                    orbits[0].append(valn)
                    break

        if (current_orbit == total_orbits) and current_orbit!=0:
            total_orbits = total_orbits + 1
            orbits.append([])
            current_orbit = total_orbits
        else:
            current_orbit = total_orbits

    for o in orbits:
        mystr=""
        for n in o:
            mystr = mystr + str(n) + "\t"
        print(mystr)
  
    mystr=""
    for o in orbits:
        for n in o:      
            mystr = mystr + str(n) + "\t"
            break
    print(mystr)
orbits()

Hobbymath's wanderings #8

Proofs of Fermat's theorem on sums of two squares.

Brahmagupta–Fibonacci identity.

Heronian triangle.

Hamiltonian path.

Cardioid.

Sophomore's dream.

Trefoil knot.

Triangle read by rows (OEIS A125605).

Proof by infinite descent.

Tossing the Prime Coin | Understanding the Riemann Hypothesis.

Divisibility of central binomial coefficients.

Magic square of squares.

Möbius function two-dimensional (pseudo)random walk vs random walks

Playing with random walks is quite fascinating, in this case instead of using a random generator I have applied the Möbius function to decide the next step in a two-dimensional walk. It is said that the behavior of the Möbius function is pseudorandom. Indeed looking at the output of the walk it looks really like a "standard" random walk. Here is the output:

And this is a random walk generated with a random seed:

I posted a question regarding "reset random walks" at MSE recently. Basically, when a reset condition is added to the random walk, sometimes there is a clear symmetry on the graph, and sometimes not, so it is not clear which reset conditions have as a result a symmetrical graph or not. Here is an example of symmetry:

Here is the link to the question: Why does symmetry happen in reset-based random walks?

And here is the Python code. Have fun!

# Mobius walk vs Random walk
# by David M.@http://hobbymaths.blogspot.jp/

def mbrw():
    from sympy import mobius
    import matplotlib.pyplot as plt
    from random import randint
   
    def calc_dir(last_dir,mob):
        if last_dir == "R":
            if mob == -1:
                return "U"
            if mob == 0:
                return "R"
            if mob == 1:
                return "D"
        elif last_dir == "D":
            if mob == -1:
                return "R"
            if mob == 0:
                return "D"
            if mob == 1:
                return "L"
        elif last_dir == "L":
            if mob == -1:
                return "D"
            if mob == 0:
                return "L"
            if mob == 1:
                return "U"
        elif last_dir == "U":
            if mob == -1:
                return "L"
            if mob == 0:
                return "U"
            if mob == 1:
                return "R"

   
    testlimit = 1000000
    lox = []
    loy = []
    last_dir = "R"
    is_first = True
    for n in range (1,testlimit):
       
        # Function used to decide the direction
        # mobius or random (uncomment for random seed)
        tmpval=mobius(n)
        #myrand = randint(0,9)
        #if myrand <=3:
        #    tmpval = -1
        #elif myrand <=6:
        #    tmpval = 0
        #else:
        #    tmpval = 1

        if is_first:
            if n==1:
                # Starting condition
                last_dir = "R"
                lox.append(0)
                loy.append(0)
                lox.append(1)
                loy.append(0)
               
            is_first = False
        else:
            # Calculation of next direction
            last_dir = calc_dir(last_dir,tmpval)
            nextx=0
            nexty=0

            # Displacement in the new direction
            jumpvalue = 1
           
            # This is a sample of symmetrical reset condition
            # uncomment to apply a reset condition
            #if is_square(n):
            #    nextx=0
            #    nexty=0
            #else:
            # Calculation of next (x,y)
            if last_dir=="R":
                nextx=lox[n-1]+jumpvalue
                nexty=loy[n-1]
            elif last_dir=="D":
                nextx=lox[n-1]
                nexty=loy[n-1]-jumpvalue
            elif last_dir=="L":
                nextx=lox[n-1]-jumpvalue
                nexty=loy[n-1]
            elif last_dir=="U":
                nextx=lox[n-1]
                nexty=loy[n-1]+jumpvalue
               
            lox.append(nextx)
            loy.append(nexty)
   
    plt.plot(lox,loy)
    plt.show()

mbrw()