$$S_{D}=\{(x,y): (x_{n+1},y_{n+1}) = (Im(\psi(x_{n} + y_n i)\cdot\sin{\frac{D}{f}}, Re(\psi(x_{n} + y_n i))\cdot\cos{\frac{D}{f}})\}$$
Where:
1. The seed is $(x_0,y_0)=(\pi/4,\pi/4)$
2. $\psi(x_{n} + y_n i)$ is the digamma function applied to a complex number $z_n=x_n+y_n i$
3. $f$ is a fraction constraint initially for this test fixed to $f = 1$.
4. $D$ is a constraint called "Depth", $D \in \Bbb N$.
5. The process is just calculating on each iteration the digamma function and then generating a complex number $(z_{n+1}=x_{n+1},y_{n+1})$, where the real part of the complex number at the iteration $n+1$ comes from the imaginary part of the complex number associated with the previous iteration $n$, and likewise the imaginary part of the complex number of the iteration $n+1$ comes from the real part of the complex number generated at the iteration $n$. In other words, it is "crossing" $Im$ and $Re$ on each iteration.
Some details of the trial and error process: Any other combination of sine, cosine or not "crossing" $Im$ and $Re$ on each iteration will not provide the patterns I will show below. Modifying the seed will provide the same results but located in other points of the plane or distorted. Modifying the fraction $f$ just generates other similar attractors, so $f$ was fixed to focus the tests.
Basically: for every fixed $D$ (Depth) and for a fixed $f=1$ value, if we run the system up to as much as we want big $n$ we arrive to a sink of a strange attractor, so from that point the set of points will not grow, and the graph of the set of points $S_n$ shows an attractor.
Running the system for $D \in [0,250]$ and up to $n = 4 \cdot 10^4$ (it can be a higher value but for the tests is a good enough number) we obtain a family of attractors that can be divided into three different types or shapes: **clusters**, **spirals** and what I call "**anomalies**".
1. Clusters are just clusters of points, a rounded shape, but the points can not escape from that space. For instance, $S_5$:
2. Spirals have the shape of different types of spirals: two arms, three arms, four, etc. For instance, $S_{45}$:
3. Anomalies are very specific shaped attractors. The limits are very well defined. For instance, $S_{22}$ and $S_{75}$:
I wanted to see all them together, they look like this, it is zoomed a little bit:
There is too much density of points when we show the big clusters together with the spirals and the anomalies, so I decided to split them into three groups and I focused on the family of spirals, and this is the graph of the spiral-shaped attractors in $Depth \in [0,250]$ (
click to zoom the image in a full screen view):
I wanted to make a simulation exercise with this. As it is said that complex dynamical systems are closer to reality and the attractors obtained are spirals and clusters I tried to see if it is possible to simulate some galaxies and star clusters with this information. So I added a third dimension using the $Depth$ value.
$$S_{Depth}=\{(x,y,d): (x_{n+1},y_{n+1},d_{n+1}) = (Im(\psi(x_{n} + y_n i)\cdot\sin{\frac{D}{f}}, Re(\psi(x_{n} + y_n i))\cdot\cos{\frac{D}{f}}, D+(\sin{n}))\}$$
The $d$ value locates each attractor $S_D$ in a different depth, and the sine is just a trick to add some "random" volume to each cloud of points of each attractor. I used a
LiDAR free viewer to visualize the $(x,y,d)$ points in the three-dimensional space (in my case
PointCloudViz Free Edition), and it looked quite promising. I decided to invert and colorize the images with GIMP2, and these are some results. They are exactly the same points as shown in the previous plots, but colorized and viewed from inside the cloud of points in a three-dimensional environment (
click to zoom the images in a full screen view):
It seems that these points can "simulate" galaxies at some rudimentary extent. It would be interesting to run a
n-body simulation with this family of attractors.
I am not sure if it is correct to talk about strange attractors in the case of the cluster and spiral shapes. The iterations really arrive to a $n$ value such that the point $(x_n,y_n)$ is a sink and the set of points does not grow from that moment, but I am not sure if it is correct to talk about strange attractor, just attractor or just dynamical system attached to a complex map. I am more sure about for instance the case $D=22$ shown above because the shape of the attractor is very clear.
I wonder if there are any other dynamical systems able to generate this kind of shapes.
update 2017/11/10: for the sake of completeness, below there is a two dimensional $n$-body gravity simulation using the first $600$ points generated by the $S_5$ set. We consider the points as a
protoplanetary disk, and the closest point to the average center of the points will be given the biggest density and thus it will turn into a central sun. I have provided random densities, radius, mass and momentum to the rest of points, so they simulate distinct types of planets by accretion. The less dense planets are whitish-bluish, the average density planets are greenish (Earth-like) and the heaviest planets are reddish. Depending on the volume, the colors are lightest on each variety, the bigger the lighter. When two planets collide, the bigger mass takes the average density of both planets and all the mass and momentum. When the planet collides with the central sun, the sun takes all the mass and the average density but it remains fixed as the center of the system. This simulation (is a looped animation) has $12500$ units of time. The Python code is a very nice work by Thanassis Tsiodras, with some customizations of my own to add random densities and colors for a little bit more realistic visualization, so credits to the author (
site here):
It seems that the clouds of points generated by the dynamical system are valid to simulate some basic behavior of protoplanetary disks. We can see that the initial accretion cloud starts to be cleaned out: some protoplanets collide and generate a bigger planet with a new averaged density and radius. Others are expelled from the protoplanetary disk due to the distinct gravity interactions. Others collide with the sun, making it bigger, and finally after some units of time it gets stable with some few planets of different density and more or less stable trajectories. Be aware that the animation is gradually zooming out to see the more external trajectories (it seems that the planets and the sun are shrinking because there is no background reference).
Update 2017/11/20: I have uploaded the code to my repl.it account (
link here). It cannot be run online because it requires saving the images of the different dynamical systems. It will work in, for instance, a personal computer using Python Anaconda. Feel free to use and modify it and enjoy! Just in case, here is a backup copy too:
# Why is this family of dynamical systems able to produce spirals and clusters of points?
# code attached to the question at Mathematics Stack Exchange: https://math.stackexchange.com/q/2506900/189215
# keys: recurrence-relations, dynamical-systems, visualization, experimental-mathematics, chaotic-systems
def sa():
import matplotlib.pyplot as plt
from math import pi
from cmath import cos , sin
from scipy.special import digamma
# Each system is iterated up to [testlimit] times
# It can be a higher value but for the tests is a good enough number to arrive to the sink points
testlimit = 40000 #1000000
# 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 [testlimit] points. If we arrive to a sink point before the [testlimit]
# 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 = []
# Depth D parameter range
rangestart = 0
rangelimit = 250
# f parameter definition
fraction=1
# Using this very code, and reviewing the intependent graphs of the result of each
# dynamical system, it is possible to split them into three different type of systems
# anomalies, big clusters and spirals, see explanation at https://math.stackexchange.com/q/2506900/189215
# you do not know beforehand (without printing) which Depth D value gives an anomalie, big cluster or spiral type
# system, so the code below was added after printing individually each system and then splitting them into the
# three possible types. So it was made by trial-error.
anomalies_list = [22,25,44,47,66,88,97,110,113,132,135,151,173,176,179,195,198,204,220]
bigclusters_list = [5,8,11,16,27,30,33,38,49,52,55,60,71,74,77,79,82,93,96,99,104,115,118,121,126,137,140,143,147,148,159,162,165,169,170,180,181,184,187,192,202,203,206,209]
# We will generate a dynamical system for each Depth D parameter from [rangestart] to [rangelimit] included.
for depth in range (rangestart,rangelimit+1):
# Uncomment to avoid plotting the anomalies
#if depth in anomalies_list:
# continue
# Uncomment to avoid plotting the big clusters
#if depth in bigclusters_list:
# continue
st=[]
# Seed of the dynamical system
st.append([pi/4,pi/4,depth])
# We will generate up to [testlimit] points of the current dynamical systes for a fixed D Depth.
# if we arrive to a sink point we do not need to continue the loop up to [testlimit] 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,testlimit):
# 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 digamma function of the previous iteration complex point [oldx+(oldy*(1j)]
c=digamma(oldx+(oldy*(1j)))
# This is the new point of the dynamical system for the current iteration
newx=(c.imag)*(sin(depth/fraction))
newy=(c.real)*(cos(depth/fraction))
# The z coordinate is just a trick to give a random depth in space to the points
# so we can simulate at some rudimentary extent a "physical" galaxy with the (x,y,z) points later
newz=depth+(sin(n)).real
# We add it to the list [st] of points (x,y,z) of the current dynamical system
st.append([newx.real,newy.real,newz])
# Once we have generated the points of the current D, f dynamical system, we have a
# set of points S_D (where D is the Depth number). E.g S_1, S_2, etc.
# see examples at https://math.stackexchange.com/q/2506900/189215
# this loop below will remove repeated points before plotting them
# to avoid memory overflow
for idx in range(0,len(st)-1):
if len(lox)==0:
lox.append(st[idx][0])
loy.append(st[idx][1])
loz.append(st[idx][2])
elif not(str(lox[len(lox)-1])==str(st[idx][0]) and (str(loy[len(loy)-1])==str(st[idx][1]))):
# if we arrive to a sink we cannot go out from it, so we do not need to add again the same point
lox.append(st[idx][0])
loy.append(st[idx][1])
loz.append(st[idx][2])
# we set black color to the background and white color to the points
ax = plt.gca()
ax.set_axis_bgcolor((0, 0, 0))
figure = plt.gcf()
# If we want to see accuracy and details, we need to use a big enough resolution like the one below
# Uncomment if you want to use the ZOOM-IN code
#figure.set_size_inches(320, 32)
#figure.set_size_inches(20, 32)
##ZOOM-IN print code BEGIN
#if depth==rangelimit:
# # 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 loop will generate some png files making a zoom-in to the plot
# # of the family of dynamical systems. You can joint them and create
# # an animation (div, anigif) with for instance the free application [VirtualDub]
# # the values and limits of the zoom were calculated by trial and error
# for zoom in range (105,140):
# dim=37/(zoom/100)
# plt.axis([-dim,dim,-dim*4/37,dim*4/37])
# plt.savefig("dynamical_system_digamma_family_zoom_in_"+str(zoom)+".png")
##ZOOM-IN print code END
# We will show the current Depth D generated dynamical system
plt.plot(lox,loy,",w")
#plt.show()
#We will save each dynamical system in a png file
plt.savefig("dynamical_system_digamma_family_depth_"+str(depth)+".png")
# then we clean the buffers and lists
# If you want to use the ZOOM-IN code bove please COMMENT the code below (BEGIN-END)
# BEGIN Clean buffers
plt.clf()
lox = []
loy = []
# END Clean buffers
sa()
