Inverse Iteration Method for Julia Sets

from Essential Mathematica for Students of Science
© James J. Kelly, 1998
Last revised: January, 2006

Introduction

The inverse iteration method (IIM) is based upon a theorem due to Julia and Fatou stating that a Julia set consists of all periodic repelling points of a nonlinear mapping function on the complex plane.  To understand how this theorem can be used as the basis for algorithms that construct Julia sets, we review some of the terminology and properties of iterated dynamical systems.

Let z -> f[z, c] represent a mapping of points in the complex plane generated by a function f with parameter c.  Forward iteration of the map beginning with seed (Subscript[z, 0]) produces a sequence of points {(Subscript[z, 0]),  (Subscript[z, 1]),  ⋯,  (Subscript[z, k])}, where \(TraditionalForm
      \`z\_k = \(f[z\_\(k - 1\), c] = 
          \(\(f\^k\)[z\_0, c] = f[\(f\^\(k - 1\)\)[z\_0, c]]\)\)\) represents the (k^th) iteration of the mapping function.  This sequence of points is known as the forward orbit of (Subscript[z, 0]).  Alternatively, the backward orbit is constructed by inverse iteration \(TraditionalForm
      \`\(f\^\(-k\)\)[z, c] = \(f\^\(-1\)\)[\(f\^\(1 - k\)\)[z, c]]\) where (f^(-1))[z, c] is the inverse of f.  However, because the inverse function for interesting systems is usually multivalued, the backward orbit has a tree structure in which the number of points grows geometrically as the iteration progresses.  This basic structure is illustrated below.  The points obtained by inverse iteration are described as pre-images of the initial point.

Any point z satisfying z == f[z, c] is a fixed point of the map.  Fixed points can be classified as
a) attracting if |(f^′)[z, 
            c]|  <  1
b) repelling if  |(f^′)[z, 
            c]|  >  1
c) neutral or indifferent if |(f^′)[z, 
            c]|  =  1
Neutral fixed points can be classified further by expressing the derivative of the mapping function as \(TraditionalForm\`f\^′ = e\^\(2  π\ i\ α\)\) where α is real.  Thus, a fixed point is described as either rationally or irrationally indifferent according to α.  Similarly, any point z satisfying z == (f^k)[z, c] and (f^j)[z, c] ≠ z for any 1 ≤ j < k is classified as a periodic point and identified as a member of a cycle of period k or a k-cycle for short.  Periodic points are either attracting, repelling, or indifferent according to value of \(TraditionalForm\`|\((f\^k)\)\^′|\).

Consider a point z near a periodic point (Subscript[z, k]).  If the cycle is repelling, application of the function (f^k) increases the distance between the test point and the periodic point, whereas for attractive cycles that distance decreases.  Conversely, inverse iteration (f^(-k)) produces a set of points that is closer to a repelling k-cycle.  Thus, according to the theorem that the Julia set consists of repelling periodic points, we expect that any seed point subjected to inverse iteration long enough will generate a set of points that approximates the Julia set.  This behavior is reminiscent of other dynamical systems, such as the Henon map, for which there is a unique strange attractor such that a point is either captured by the set or escapes to infinity upon iteration.

In this notebook we apply the IIM to the quadratic map f[z, c] = (z^2) + c and leave other maps to the exercises.   For this map the inverse is double-valued, (f^(-1))[z, 
          c] = ±(Sqrt[z - c]).  Given a seed on (Subscript[z, 0]) ∈ J[f, c], the first inverse iteration yields 2 points, the next 4 points, and the (k^th) iteration (2^k) more points approximating J[f, c].  Thus, the deterministic version of the IIM algorithm uses Nest to produce list of (2^k) points representing the (k^th) inverse iteration beginning from an arbitrary seed.  Alternatively, one can construct a set of points by selecting one of the two roots at random.  Some number of initial iterations, depending upon c and the seed chosen, must be discarded to allow transients to decay, but all subsequent pre-images are retained.  For a given number of points in the final set, the random IIM reaches much deeper into the iteration tree than the deterministic algorithm and, hence, usually converges more rapidly.

References

The Science of Fractal Images, edited by H.-O. Peitgen and D. Saupe, (Springer-Verlag, Berlin, 1988).

H.-O. Peitgen and P.H. Richter, The Beauty of Fractals, (Springer-Verlag, Berlin, 1986).

Initialization

ClearAll["Global`*"];
Off[General::spell,  General::spell1];
$TextStyle = {FontFamily -> "Times", 
          FontSize -> 12};

Needs["Graphics`"]

ComplexToArgand = {Complex[x_, y_] -> {x, y}, 
          x_Real -> {x, 0}, x_Integer -> {x, 0}};

Implementation

Deterministic IIM algorithm

The deterministic IIM algorithm replaces each point by two pre-images and iterates the procedure n times to produce (2^n) points approximating the Julia set.  The simplest implementation uses Nest to perform the iteration and Flatten to produce a simple list of complex numbers.  To ensure that all results are numerical instead of symbolic, we convert arguments to numerical form before iteration, which saves time and storage.  The final results are then converted from complex numbers to Argand coordinates for plotting.

\(\(\(JuliaDIIM::usage\  = \ \*""JuliaDIIM[c,z0,n] performs n inverse iterations of the quadratic map for parameter c beginning with seed z0, producing (2^n) points approximating the Julia set."";\)\(
\)
    \)\)
JuliaDIIM[c_, z0_, n_] := 
      Module[{cc = N[c], 
          zz = N[z0]}, (Nest[{(Sqrt[# - cc]), -(Sqrt[# - cc])} &, zz, n] // 
              Flatten) /. ComplexToArgand]

Note that this implementation retains only the final level of the iteration tree.

Random IIM algorithm

The random IIM algorithm chooses one of the two pre-images at random, effectively performing a random walk on the inverse iteration tree.  All hits on this tree are retained using NestList.

\(\(\(JuliaRIIM::usage\  = \ "JuliaRIIM[c,z0,n] samples n points of the backward orbit of seed z0 for the quadratic map with parameter c, choosing the root randomly.";\)\(
\)
    \)\)
JuliaRIIM[c_, z0_, n_] := 
      Module[{cc = N[c], zz = N[z0]}, 
        NestList[(1 - 2  Random[Integer]) (Sqrt[# - cc]) &, zz, n] /. 
          ComplexToArgand]

Note that a number of the initial iterations, depending upon the seed, should be discarded before plotting in order to suppress transients before the inverse orbit has settled on the attractor.

Histograming frequency distributions

The IIM algorithms often visit various portions of a Julia set with vastly different frequencies, hitting some regions often and others rarely.  Consequently, for some sets it can take a very large number of iterations before an adequate representation of the set is achieved.  This behavior can be investigated using the following function which produces a two-dimensional table of hit frequencies on a uniform grid.  The resulting two-dimensional histogram can be displayed using ListPlot3D, for example.

histogram[
        data_, {xmin_, xmax_, nx_Integer}, {ymin_, ymax_, ny_Integer}] := 
      Module[{histo, i, j, k, dx, dy}, dx = (xmax - xmin)/nx; 
        dy = (ymax - ymin)/ny; histo = Table[0, {nx}, {ny}]; 
        Do[{i, j} =  
            Round[(data[[k]] - {xmin, 
                      ymin})/{dx, dy}]; i = Min[nx, Max[i, 1]]; 
          j = Min[ny, Max[j, 1]];  ++histo[[i, 
              j]], {k, 1, Length[data]}]; 
        histo]

Examples

You may find that in the interest of economizing file size, this notebook may have been provided without output, but regeneration of the output does not take long on modern microcomputers.

Sampling properties

We begin with a connected Julia set belonging to the main bud of the Mandelbrot set.

JuliaData = JuliaDIIM[-1, 0.2, 15];

ListPlot[JuliaData, Axes -> False, 
        AspectRatio -> Automatic, 
        Prolog -> AbsolutePointSize[1]];

JuliaData = 
        Drop[Evaluate[JuliaRIIM[-1, 5  I, 100000]], 100];

ListPlot[JuliaData, Axes -> False, 
        AspectRatio -> Automatic, 
        Prolog -> AbsolutePointSize[1]];

Both algorithms produce very similar images and the choice of seed does not seem to matter (although there are special cases which should be avoided, such as (Subscript[z, 0]) = 0).

Upon closer examination one finds that some regions are rather sparsely populated.

ListPlot[JuliaData, Axes -> False, 
        AspectRatio -> Automatic, Prolog -> AbsolutePointSize[1], 
        PlotRange -> {{-0.12, 0.12}, {0.5, 0.8}}];

This nonuniform sampling is conveniently displayed in the form of a two-dimensional histogram of hit frequencies.

histo = 
        histogram[
          JuliaData, {-1.6, 1.6, 161}, {-0.8, 0.8, 81}];

ListPlot3D[histo, Boxed -> False, 
        Axes -> False];

The IIM does not work well for the case below — comparison with the results of the escape time algorithm shows that IIM misses practically of the interesting structure within the outline shown unless a prohibitively large number of points is employed; see Julia.nb for comparison.  Nevertheless, for dendritic sets IIM is much more efficient than the escape time algorithm.  More sophisticated algorithms are available which are more successful than either for other sets.

JuliaData = 
        Drop[Evaluate[JuliaRIIM[-0.745429 + 0.113008  I, I, 100000]], 
          100];
ListPlot[JuliaData, Axes -> False, 
        AspectRatio -> Automatic, 
        Prolog -> AbsolutePointSize[1]];

Fancier algorithms keep track of the hit frequencies and stop accumulating points in regions visited too often; in that way the proportion of samples in less popular regions is enhanced.  However, we will not pursue such algorithms here.

Dendrites

IIM algorithms are usually better for dendritic sets than escape time or boundary scanning algorithms, in part because they do not waste a lot of time examining large areas not relevant to the desired Julia set.

JuliaData = Drop[Evaluate[JuliaRIIM[I, 1, 100000]], 100];
ListPlot[JuliaData, Axes -> False, 
        AspectRatio -> Automatic, 
        Prolog -> AbsolutePointSize[1]];

It is interesting to examine the convergence of the deterministic algorithm given a distant seed.  If we begin at a distant point and perform too few iterations, the set appears to enclose significant empty space

JuliaData = JuliaDIIM[I, 10, 10];
ListPlot[JuliaData, Axes -> False, 
        AspectRatio -> Automatic, 
        Prolog -> AbsolutePointSize[1]];

but upon further iteration the enclosed space collapses and reveals the filamentary nature of this set.

JuliaData = JuliaDIIM[I, 10, 15];
ListPlot[JuliaData, Axes -> False, 
        AspectRatio -> Automatic, 
        Prolog -> AbsolutePointSize[1]];

JuliaData = JuliaDIIM[I, 10, 17];
ListPlot[JuliaData, Axes -> False, 
        AspectRatio -> Automatic, 
        Prolog -> AbsolutePointSize[1]];

Notice that a different starting condition produces more rapid convergence.

JuliaData = JuliaDIIM[I, 0, 10];
ListPlot[JuliaData, Axes -> False, 
        AspectRatio -> Automatic, 
        Prolog -> AbsolutePointSize[1]];

and that the random IIM yields good convergence without misleading enclosed space even for a distant seed.  Below we choose the number of points to be about (2^10), plus some discarded initial iterations, for comparison with the 10 iteration DIIM result at the beginning of this section.

JuliaData = Drop[Evaluate[JuliaRIIM[I, 10, 1200]], 100];
ListPlot[JuliaData, Axes -> False, 
        AspectRatio -> Automatic, 
        Prolog -> AbsolutePointSize[1]];

The RIIM algorithm is more robust because it penetrates deeper into the inverse iteration tree for the same total number of points with points sprinkled over all levels, whereas the DIIM produces all points at a single shallower level.