Continuous iteration of fractals.

(Extra picture for search engine related reasons)

Like many others, I was intrigued when I first saw pictures of Julia sets and Mandelbrod sets, and I had lots of fun generating them on a computer.

But does this have anything to do with physics? Is there a dynamic system whose behavior is modeled by the startlingly simple equation that so surprisingly 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:

  z->f(z)->f2(z)->f3(z)-> ... ->fN(z) = z^2 + c

For large N, we would start to get something like a dynamic system in continuous time.
Taking the Nth iterative root of a function is apparently a difficult subject. I’ll just list a few cases where the answer is simple.

The map
   z-> z + c
has as Nth root
    z-> z+c/N

The map
   z-> az
has as Nth root
    z-> (a^(1/N)) z

Note that already, we come across multi-valuedness: there are N Nth roots of a complex number. We will come back to multi-valuedness later.

Mobius transformations :
  z-> (az+b)/(cz+d)
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 structure as functions under composition, is an interesting approach. An online article 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 matrix form:


     (z^0)    ( 0    0   0   0 ...........)  (z^0)  
     (z^1)    (a0   a1  a2  a3............)  (z^1)  
d/dt (z^2) =  ( 0  2a0 2a1 2a2 2a3 .......)  (z^2)   
     (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 functions 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 one, we need the higher columns and rows to get progressively less important. But in our case, they don’t.

At this point, I found the web site of Markus Mueller . Apparently, he had succeeded in making continuous iterations of the real case (Verhulst dynamics). He has some nice pictures of it on his web page, 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 iteration has an attractive fixed point (z0). Near this fixed point, it becomes increasingly well approximated by linear dynamics:

   (z-z0) -> lambda * (z-z0)

So near this point, we can find continuous iterates. The nice idea that Markus Mueller had was to then simply forward iterate back to the original point, using “ordinary” iteration.

Because Markus Mueller uses real numbers, the multi-valuedness of the square roots and the Nth roots of lambda are not so worrying. Each choice of either the positive square root or the negative square root leads to a different fixed point. We choose a combination with a small real-valued Nth root of lambda, since that’s the only one that will give smooth real-valued continuous iterates.

This method appears to work quite well numerically.

However, the time evolution generated in this way apparently cannot be that of a dynamical system in 1 variable. This can be seen by the fact that for example it crosses the line z=0 may times in time, but each time follows a different trajectory afterwards: The value of z(t) itself does not define a state from which z(t+dt) is uniquely determined.

I decided to try the method on complex numbers instead of just reals. After all, it is the complex numbers which give the nice Julia pictures. Again, we run into a multi-valuedness problem. To produce some pictures, I initially just ignored the multi-value problem, by simply choosing roots. Again, the method works numerically, and soon I was having Lissajous-like curves on my screen, that winded through the complex plane. But unfortunately, there were usually self intersections. These intersections are bad, because they 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 linked to the multi-valuedness: If all maps and their inverses  the story were single-valued, then we should not get self intersections.

There is a way to get rid of multi-valuedness: Use a Riemann surface! In 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 u=(R,theta). These are just polar coordinates of the complex plane, except that theta 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 2pi.

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 plane, adding a constant is just a translation. If we construct the Riemann surface of unwrapped-C, we make a number of copies of C, cut them open along a line from infinity to the origin, and the glue them together like multi-story spiraling parking garage. On this surface, we can do a translation by a constant (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 it by c, we see that the small region gets sliced in 2: each half ending in a different floor  of the parking garage.
We will just accept the fact there is a discontinuity in this map, and define z->z+c as a translation in the local copy of C within unwrapped-C.

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 on 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, and 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 will evolve away from it after an infinitely small perturbation. In our coordinates, the dynamics looks like an explosion: The system goes away from z=1 at ever increasing rate. That is, using logarithmic scale. If we used linear scale, 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 region in the plane. One thing that was surprising to me, is that the border between the interior and exterior is not a separatrix. This might have been a problem, because the border is an infinitely long line (like the coast of Britain). 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 a fixed point! The clever way that the system fulfills both, is by becoming an oscillation of infinite amplitude, but sampled at a certain sample rate, can appear to be always finite at each snapshot!! In other words, an infinite oscillation that looks finite when viewed through a stroboscope. Check that rectangles on lines that cross the border are either all outside or all inside the interior.

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 dynamics, even some regions that appear to stay smooth for all time. Not only the lines 0-(-c) are discontinuous, also image of those lines under z->z^2+c. So 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 discontinuities 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 extraordinary properly is "stroboscopy", infinite oscillation that appears finite when 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.