Continuous iteration of fractals.
picture for search engine related reasons)
Like many others, I was intrigued when I first saw pictures of Julia
and Mandelbrod sets, and I had lots of fun generating them on a
But does this have anything to do with physics? Is there a dynamic
whose behavior is modeled by the startlingly simple equation that so
generates all those pretty pictures:
z -> z^2 + c
To get a continuous time evolution, we need the “Nth
iterative root” of this map, a map f(z), such that:
->fN(z) = z^2 + c
For large N, we would start to get something like a dynamic system in
Taking the Nth iterative root of a function is apparently a difficult
I’ll just list a few cases where the answer is simple.
z-> z + c
has as Nth root
has as Nth root
z-> (a^(1/N)) z
Note that already, we come across multi-valuedness: there are N Nth
of a complex number. We will come back to multi-valuedness later.
form a group, and they are a representation of the group of 2X2 matrices
( a b )
( c d )
So taking the Nth root of this matrix, will produce coefficients for an
Nth root Mobius transformation.
Actually, finding a group of matrices which have the same group
as functions under composition, is an interesting approach. An online
on this is by
Aldrovandi and Freitas
Consider a non-linear dynamical system with 1 variable:
dz/dt = f(z)
Expand f(z) into a Taylor polynomial:
f(z) = a0z^0 + a1z^1 + a2z^2 + ....
Using some fairly easy differentiation rules, we can write this in a
(z^0) ( 0
0 0 0 ...........)
a1 a2 a3............) (z^1)
d/dt (z^2) = ( 0 2a0 2a1 2a2 2a3 .......)
(z^3) ( 0 0
3a0 3a1 3a2 3a3 ...) (z^3)
In a slightly modified form this matrix is called the Bell matrix of
the function f.
It turns out that Bell matrices form a group, and that compositions of
correspond to multiplications of their Bell matrices and vice versa.
So all we need to do is find the Nth root of the Bell matrix, and we
will have the required Nth iterative root!
The problem is, that to truncate the infinite Bell matrix into a finite
we need the higher columns and rows to get progressively less
But in our case, they don’t.
At this point, I found the web site of Markus
. Apparently, he had succeeded in making continuous iterations of the
case (Verhulst dynamics). He has some nice pictures of it on his web
and he was kind enough to explain his method to me.
First, he reverse-iterates a point a number of times by using the
inverse iteration: z->sqrt(z-c). It turns out the reverse
an attractive fixed point (z0). Near this fixed point, it becomes
well approximated by linear dynamics:
(z-z0) -> lambda * (z-z0)
So near this point, we can find continuous iterates. The nice idea that
Mueller had was to then simply forward iterate back to the original
using “ordinary” iteration.
Because Markus Mueller uses real numbers, the multi-valuedness of the
roots and the Nth roots of lambda are not so worrying. Each choice of
the positive square root or the negative square root leads to a
fixed point. We choose a combination with a small real-valued Nth root
lambda, since that’s the only one that will give smooth
This method appears to work quite well numerically.
However, the time evolution generated in this way apparently cannot be
of a dynamical system in 1 variable. This can be seen by the fact that
example it crosses the line z=0 may times in time, but each time
a different trajectory afterwards: The value of z(t) itself does not
a state from which z(t+dt) is uniquely determined.
I decided to try the method on complex numbers instead of just reals.
all, it is the complex numbers which give the nice Julia pictures.
we run into a multi-valuedness problem. To produce some pictures, I
just ignored the multi-value problem, by simply choosing roots. Again,
method works numerically, and soon I was having Lissajous-like curves
my screen, that winded through the complex plane. But unfortunately,
were usually self intersections. These intersections are bad, because
means that we do not have a dynamical system uniquely defined by z(t).
After some thought, I realized that he self-intersections were in fact
to the multi-valuedness: If all maps and their inverses the
single-valued, then we should not get self intersections.
There is a way to get rid of multi-valuedness: Use a Riemann surface!
our case, we need a Riemann surface that I will call
“unwrapped-C”. It is
basically the Riemann surface of the complex logarithm.
In unwrapped-C, each number (u) can be written as 2 real numbers
These are just polar coordinates of the complex plane, except that
is no longer limited to a 2pi interval.
To multiply in unwrapped-C, we use
u1*u2 = (R1*R2,theta1+theta2)
This is the same as in ordinary C, but we retain the extra multiples of
There is now a unique Nth root, or pth power:
u1^p = (R1^p, p*theta)
The tricky part is now to add a constant (c). On the ordinary complex
adding a constant is just a translation. If we construct the Riemann
of unwrapped-C, we make a number of copies of C, cut them open along a
from infinity to the origin, and the glue them together like
spiraling parking garage. On this surface, we can do a translation by a
(c) just as in C. But look at the line extending from –c to
the origin. If
we take a small region that contains part of this line, and translate
by c, we see that the small region gets sliced in 2: each half ending
a different floor of the parking garage.
We will just accept the fact there is a discontinuity in this map, and
z->z+c as a translation in the local copy of C within
OK, so we now have a map z->z^2+c whose inverse is single
valued, at the cost of a discontinuity. So lets make pictures!
The pictures I made use a coordinate system in which the horizontal
axis is theta, and the vertical axis is log(R). If we draw a Julia set
it, we will get a periodic wall paper of deformed ordinary Julia sets.
Since continuous iterations will involve back-iterating to the fixed
point, I thought we might as well start near the fixed point. The fixed
points are small red dots. Around the fixed point, I make a small
circle of starting M points.
I then iterate these points N times, using the iterative Nth root of a
“normal” iteration. These N points are joined in a
line. The M lines of N points are then normally iterated. On each line,
I also put 1 small rectangle, so that you can see what a
“normal” iteration does: it jumps from rectangle to
rectangle along a curve.
The blue lines can be interpreted as time evolutions of the stateof the
system in continuous time.
I also insert red lines going from –c to the origin. If a
curve hits this line, it will “jump”.
So what do we get?
The first example is c=0, so z-> z^2.
(Click to enlarge)
This case can also be solved exactly. There is a fixed point at z=1,
a fixed point at z=0. The point z=0 is grotesquely deformed at
-infinity (of the vertical scale),
because of the use of logarithm. The point z=1 is unstable, the state
evolve away from it after an infinitely small perturbation. In our
the dynamics looks like an explosion: The system goes away from z=1 at
increasing rate. That is, using logarithmic scale. If we used linear
we would see that the lines in the black part of the Julia picture all
move to z=0.
Next we try to slightly move c. We use c=-0.2.
(Click to enlarge)
We again get a seemingly well behaved picture, at least for a certain
in the plane. One thing that was surprising to me, is that the border
the interior and exterior is not a separatrix. This might have been a
because the border is an infinitely long line (like the coast of
But we get orbits that cross the border, and points in the exterior
are supposed to go to infinity, while those on the interior attract to
point! The clever way that the system fulfills both, is by becoming an
of infinite amplitude, but sampled at a certain sample rate, can appear
be always finite at each snapshot!! In other words, an infinite
that looks finite when viewed through a stroboscope. Check that
on lines that cross the border are either all outside or all inside the
Next is an example with c = -1.
(Click to enlarge)
This time, we are getting some trouble from the red lines
attached to copies
of –c. But there are still quite large regions with smooth
some regions that appear to stay smooth for all time. Not only the
0-(-c) are discontinuous, also image of those lines under
the discontinuities get worse and worse, although certain regions
appear to be safe from them.
Finally an example with c = 0.012+i0.66
(Click to enlarge)
This value of c lies outside the Mandelbrod set. We therefore know
that images of the point c will end up in infinity. So images of the
across our red lines will presumably also end up at infinity. It seems
that in this case, there may be no safe place from discontinuities.
It seems that the dynamical systems we find, are pretty weird. An
properly is "stroboscopy", infinite oscillation that appears finite
viewed through a stroboscope.
It is possible that a better picture of the dynamics can be found by
finding a more suitable Riemann surface to plot it on.