Wednesday, April 18, 2018

Studying discrete-time dynamical systems (VII): some updates and Python code

A kind reader of Hobbymaths asked me to publish the code to visualize the images that I showed at my former post (here). So here it is, it is a bit long, but I think the comments in the code are self-explanatory. You will need Python Anaconda in order to run properly the code and generate the images. Basically you will need to save the code into a text file, for instance "dynamical_system_test.py" and then run in command line (my testing environment is Windows 10 + Python) "python dynamical_system_test.py". Currently I am preparing a paper about this topic, but it will take some months to be finished, maybe there will be something ready for the end of $2018$.

The code as it is written will generate $100$ images, which are the visualizations of the evolution of the control parameter time $t$ of the following families of dynamical systems:

$$S_{D}=\{(x,y, f(z)): (x_{n+1},y_{n+1}) = (Im(f(x_{n} + y_n  i)\cdot\sin{D}, Re(f(x_{n} + y_n i))\cdot\cos{D})\}$$

Where $f(z)=\frac{1}{(z+1)^t}$, and initially $t \in \Bbb [1,2]$ (it can be greater than $2$ but for the purpose of the code I have set $2$ as the maximum value to be calculated). It will generate the attractors for $t=1.01, t=1.02, \cdots , t=2.00$. For each $t$, the control parameters $D$ and $n$ will run through $D \in [0,2\pi]$ in $8 \cdot 10^3$ steps, and $n \in [0,8 \cdot 10^3]$, and the images are showing the region $z=[+/-5]+[+/-5]i$, for the initial condition $(x_0=0,y_0=0)$. 

Please be aware that generating the images takes time! each one of them takes around $10$ minutes in a good computer. And also you need to re-dimension the width of the images with a program like e.g. Gimp2 (the explanation is written in the code). So it will generate (between others) the following images one by one sequentially:

$f(z)=\frac{1}{(z+1)}$ (this will be the first one it generates):


$f(z)=\frac{1}{(z+1)^{1.1}}$:


$f(z)=\frac{1}{(z+1)^{1.2}}$:


$\cdots$

$f(z)=\frac{1}{(z+1)^{2}}$ (this image is the last one that will generate):


And if you use an application to join the images you will obtain the following video (it took me some days to render it):


So here is the code, enjoy it and, if you find new interesting variations of the model please let me know (please keep the reference to this post in the comments so other people can find the original information here, thanks!):

# Studying discrete-time dynamical systems (VI): mathematical models
# for symmetry and asymmetry and some biological coincidences
# keys: recurrence-relations, dynamical-systems, visualization,
# experimental-mathematics, chaotic-systems
# Done by Hobbymaths:
# http://hobbymaths.blogspot.jp/2018/01/studying-discrete-time-dynamical.html

def nds():
    import matplotlib.pyplot as plt
    from math import pi, e
    from random import randint
    from cmath import sqrt, sin , cos, tan
    from scipy.special import digamma
   
   
    # Each system is iterated up to [ntestlimit] times
    # value of n: It can be a higher value but for the
    # tests is a good enough number to arrive to the sink points
    ntestlimit = 8000
   
    # Complete list of points of each system, the first point
    # of the list is the (x,y,z) seed
    # each dynamical system will have up to [ntestlimit] points.
    # If we arrive to a sink point before the [ntestlimit]
    # iteration there is not need to continue, so the quantity
    # of elements will be smaller.
    st=[]
   
    # Splitted lists of x,y and z coordinates of the points obtained
    # on each iteration of the system stored at the [st] list
    lox = []
    loy = []
    loz = []
   
    # D parameter range
    MinD = 0
    MaxD = 8000
   
    for expo in range (100,201):
        t = expo/100
        print("Exp " + str(t))
       
        lox = []
        loy = []
        loz = []
        plt.clf()
   
        # We will generate a dynamical system for each D parameter
        # from [MinD] to [MaxD] included.
        for f in range (MinD,MaxD):
            D=(f/MaxD)*(2*pi)
           
            st=[]
           
            #Initial condition of the dynamical system
            icx=0
            icy=0
           
            st.append([icx,icy,0])
           
            # We will generate up to [ntestlimit] points of the current
            # dynamical systes for a fixed D.
            # if we arrive to a sink point we do not need to continue the
            # loop up to [ntestlimit] because the
            # point will be repited again (so we do not need to plot more
            # than once the sink point)
            for n in range (1,ntestlimit):
                # We read the former the previous iteration point
                # (oldx,oldy,oldz)
                oldx=st[len(st)-1][0]
                oldy=st[len(st)-1][1]
                oldz=st[len(st)-1][2]
               
                # Then we calculate the f(z) function for the current
                # (x_n-1,y_n-1) values.
                c=1/((oldx+(oldy*(1j)))+1)**(t)
               
                # calculate new (x_n,y_n)
                newx=(c.imag)*sin(D)
                newy=(c.real)*cos(D)
               
                # The z coordinate is just a trick to give a depth in
                # space to the points according to D
                st.append([newx.real,newy.real,D])
               
            # Once we have generated the points of the current D
            # dynamical system, we have a set of points S_D
            # (where D is the D number). E.g S_1, S_2, etc.
           
            # This loop below will detect where the cycle begins BEGIN
            tlox=[]
            tloy=[]
            tloz=[]
           
            there_is_cycle=False
           
            for idx in range(0,len(st)-1):
                if not(abs(st[idx][0])<5 and abs(st[idx][1])<5):
                    continue
                if len(tlox)==0:
                    tlox.append(st[idx][0])
                    tloy.append(st[idx][1])
                    tloz.append(st[idx][2])
                else:
                    tmpl=st[idx+1:]
                    if st[idx] in tmpl:
                        tmpval=tmpl.index(st[idx])+1
                        there_is_cycle=True
                        break
                       
                    tlox.append(st[idx][0])
                    tloy.append(st[idx][1])
                    tloz.append(st[idx][2])                       
                       
            if there_is_cycle:
                lox.extend(tlox)
                loy.extend(tloy)
                loz.extend(tloz)
                pass
            else:
                lox.extend(tlox)
                loy.extend(tloy)
                loz.extend(tloz)
                pass
            # This loop below will detect where the cycle begins END
           
            # we set black color to the background and white color to the points
            ax = plt.gca()
            #ax.set_axis_bgcolor((0, 0, 0)) # old Python code
            ax.set_facecolor((0, 0, 0))     # new Python code (change depending of Python version)
            figure = plt.gcf()
           
            print(str(f)+"\r", end='')
           
            if f==MaxD-1:
           
                #image resolution 32000x3200 pixels
                #you will need to redimension width to 3200 pixels to see the image
                #e.g. with the application Gimp2
                figure.set_size_inches(320, 32)
               
                # If we have generated all the family of dynamical systems,
                # then we will print them all together
                plt.plot(lox,loy,",w")
               
                # Uncomment if you just want to diplay them
                #plt.show()
               
                # This will save the image associated with the dynamical system in a file
                # You can create a video using VirtualDub.exe
                plt.savefig("dynamical_system_inv_1_exp_"+str(expo)+"_D"+str(MaxD)+"_n"+str(ntestlimit)+"_icx_icy_"+str(icx)+"_"+str(icy)+"_Acot5_5.png")
               
nds()



Update and caveat: I doubted about adding this to the post, because it is quite pareidolic and can be misunderstood, but just for fun! It is the following family of dynamical systems:

$$S_{D}=\{(x,y, f(z)): (x_{n+1},y_{n+1}) = (Im(f(x_{n} + y_n  i)\cdot |\sin{D}|, Re(f(x_{n} + y_n i)) \cdot |\cos{D}|)\}$$

Where $f(z)=\frac{1}{(z+1)^2}$ (so $t=2$ is fixed). As usual, the control parameters $D$ and $n$ will run through $D \in [0,2\pi]$ in $8 \cdot 10^3$ steps, and $n \in [0,8 \cdot 10^3]$, and the images are showing the region $z=[+/-5]+[+/-5]i$, for the initial condition $(x_0=0,y_0=0)$.  The only difference from the inital model is the absolute value applied to the sine and cosine functions, just that.

So this is the generated image:

Now if we rotate left $90$ degrees, and zoom in the center of the image, we have a peculiar pareidolic version of the Kubrick's movie "2001: A Space Odyssey" starchild (well, sort of, at least in my case when I saw it I remembered the reference to the movie), so for your viewing pleasure:


4 comments:

  1. THANK YOU TREMENDOUSLY FOR THIS POST!! I am interested in mathematical models of beautiful geometries. My undergraduate is in neurobiology so I love neural networks and biological simulations from computer science (flocking, swarming, cellular automata, L-systems etc)... I went to graduate school for architecture so I have a fascination for interesting mathematical geometries that may become sculptural (I keep a tumblr where I collect inspiration, including links to your posts: http://ecoarch.tumblr.com/archive) and currently I'm studying computer science. I would love to know more about dynamical systems, chaos theory, lorenz systems especially if they generate beautiful visuals, both technical papers and general information! Thank you so much!

    ReplyDelete
    Replies
    1. It seems we share same interest, btw I am a BCompSc too. Ok! then I will suggest you to read these docs:

      (1) http://www.scholarpedia.org/article/Dynamical_system

      (2) https://hypertextbook.com/chaos/strange/

      (3) http://www.algosome.com/articles/strange-attractor-identification.html

      (4) http://zitogiuseppe.com/strange.html

      (5) http://paulbourke.net/fractals/lyapunov/

      And specially you need to read about Bifurcation theory, it is very important to characterize any system you find:

      (6) https://en.wikipedia.org/wiki/Bifurcation_theory

      The closest "organic-like" attractors to the ones I found are the classical Gumowski-Mira Attractors:

      (7) http://petervandernoord.nl/blog/?p=289

      Read about Hyperchaos:

      (8) http://www.scholarpedia.org/article/Hyperchaos

      And this nice paper about attractors:

      (9) http://www.math.harvard.edu/~knill/teaching/mathe320_2014/blog/RuelleIntelligencer.pdf

      (10) I had a look to your tumblr, thank you for adding my images and your words. If you like architecture, then I guess you know Michael Hansmeyer's organic columns:

      https://youtu.be/pACNKOLZizE

      By the way... I was doubting to add it, because it could be misunderstood, but as a present I have updated the current post with an image that I have generated with my models, I call it the "starchild" like in Kubrick's movie "2001: A Space Odyssey"... it is highly "pareidolic" once you see it. Is a curiosity, but I guess you will enjoy it.

      Delete
  2. This comment has been removed by the author.

    ReplyDelete
  3. Hi I have been following your work especially with prime numbers. I myself am working heavily with prime numbers and patterns within them. I was wondering if you would be interested in collaborating and if so what would be the best way to reach you to talk more.

    ReplyDelete