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:

So this is the generated image: