Complex dynamics

Exploration with Mathematica

MandelbrotSetPlot[{-2-0.07I,-1.3+0.07I}, Frame->False, 
ImageResolution->2100,ImageSize->700, MaxIterations->250,
ColorFunction->"TemperatureMap"]

This Mathematica notebook describes how to explore complex dynamics with Mathematica. We assume that the reader knows the basics of iteration on the complex plane and how that leads to Julia sets as well as parameter spaces.

Built in tools

As of V10, Mathematica has several built in tools to generate images of Julia sets and the Mandelbrot set. They are ridiculously easy to use.

MandelbrotSetPlot

Here is the very simplest way to call MandelbrotSetPlot.

MandelbrotSetPlot[]

You can zoom in by specifying opposing corners. You can customize it further with options. Here are all the options of MandelbrotSetPlot that are not already options of Graphics.

Complement[Options[MandelbrotSetPlot],Options[Graphics]]

Let’s use it.

MandelbrotSetPlot[{-0.654+0.482I,-0.645+0.474I}, 
MaxIterations->400, ColorFunction->"StarryNightColors",
ImageSize->600, ImageResolution->1200]

Of course, it’s hard to find that image without an interactive tool to choose a zoom. We’ll see such a thing in a bit.

A more important issue with MandelbrotSetPlot is that it is very specific.

JuliaSetPlot

JuliaSetPlot is awesome. JuliaSetPlot[c] automagically generates the Julia set of (Subscript[f, c])(z)=(z^2)+c. Here’s a few:

Grid[Map[JuliaSetPlot,{{-1,-0.77+0.22I},{I,0}},{2}]]

There’s quite a little bit in that line of code - pretty cool, really. Like MandelbrotSetPlot, there are quite a few options. Sometimes, we might need to increase the iteration count just because of the nature of the Julia set.

JuliaSetPlot[0.251, MaxIterations->1000]

You might try that one without the option. Even better, you can compute the Julia set of just about any polynomial. You just specify the polynomial and its variable. Here’s a groovy cubic.

JuliaSetPlot[z^3-z^2+1,z]

There are some things that JuliaSetPlot does that we haven’t talked about much. It can compute the Julia set of a rational function.

z^3+(1/10)/z^3//Together

JuliaSetPlot[(10z^6+1)/(15z^3),z]

JuliaSetPlot also implements an algorithm called inverse iteration.

JuliaSetPlot[z^3-z-1,z,ColorFunction->None]

Complex iteration and Julia sets

The squaring function

We now examine the iteration of f(z)=(z^2) on the complex plane.  It’s not too hard to see that, if \[LeftBracketingBar](Subscript[z, 0])\[RightBracketingBar]<1, then the orbit of (Subscript[z, 0]) under the squaring function converges to zero.  If \[LeftBracketingBar](Subscript[z, 0])\[RightBracketingBar]>1, then the orbit of (Subscript[z, 0]) diverges to ∞.  The two sets, {z:\[LeftBracketingBar]z\[RightBracketingBar]<1} (the interior of the unit circle and {z:\[LeftBracketingBar]z\[RightBracketingBar]>1} (the exterior of the unit circle) form a stable set under iteration of f, also called the Fatou set.  The dynamics right on the unit circle are very unstable and highly complicated.  The unit circle is called the Julia set for this function.

Here’s a little interactive tool that let’s you play with the initial seed for the squaring function.

circlePic = Graphics[{
{Thick,Circle[{0,0},1]}
}];
Clear[f];
f[z_] = z^2;
Manipulate[
Module[{orbit},
orbit = NestWhileList[f,z0.{1,I},0.001<Abs[#]<10&,1,100];
Show[{circlePic,Graphics[{
{Line[{Re[#],Im[#]}& /@ orbit]}
}]},Axes->True, PlotRange -> 2]],
{{z0,{0.95,0.05}},{-2,-2},{2,2},Locator}]

If we start on the unit circle, we stay on the unit circle.  For most such points, the orbit bounces around aimlessly.

pts = {Re[#],Im[#]}& /@ NestList[f,Exp[I],1000];
Graphics[Point[pts], Axes->True]

Yet, there’s an underlying structure!

Graphics[{Opacity[0.2],Line/@Partition[pts,2,1]}, Axes->True]

The fundamental dichotomy for quadratic polynomials

Let’s generalize this a bit to quadratic polynomials of the form (Subscript[f, c])(z)=(z^2)+c.  In fact, much of the elementary exposition of complex dynamics is based on this family.

The dichotomy

Theorem (The quadratic escape criterion):
Suppose that \[LeftBracketingBar]z\[RightBracketingBar]≥\[LeftBracketingBar]c\[RightBracketingBar] and \[LeftBracketingBar]z\[RightBracketingBar]>2. Then (Subscript[f, c]^n)(z)->∞ as n->∞.

Proof: Since \[LeftBracketingBar]z\[RightBracketingBar]>2, we can choose λ∈ℝ such that 1<λ<\[LeftBracketingBar]z\[RightBracketingBar]-1. Then,

|(Subscript[f, c])(z)|≥((|z|)^2)-|c|≥((|z|)^2)-|z|=|z|(|z|-1)>λ|z|.

Applying this recursively, we see that \[LeftBracketingBar](Subscript[f, c]^n)(z)\[RightBracketingBar]>(λ^n)\[LeftBracketingBar]z\[RightBracketingBar]->∞ as n->∞.\[Square]

Corollary:
If we iterate (Subscript[f, c]) from a seed (Subscript[z, 0]), then one of the following must happen:

The orbit of (Subscript[z, 0]) diverges to infinity or

The orbit of (Subscript[z, 0]) stays bounded

It’s worth pointing out that this exact same dichotomy holds for the iteration of any complex polynomial.  We can make a slightly stronger statement for the quadratic family (Subscript[f, c]), however, as we know we have a uniform bound of 2.

The Julia set

Definition:
Given a polynomial f, the set of all z such that the orbit of z stays bounded under iteration of f is called the filled Julia set of f and denoted by (Subscript[K, f]).  The boundary of (Subscript[K, f]) is called the Julia set of f.

For example, the filled Julia set for the squaring function is closed unit disk.  The Julia set for the squaring function is the unit circle.

An algorithm

The fundamental dichotomy provides an algorithm to draw a filled Julia set for a quadratic function (Subscript[f, c])(z)=(z^2)+c.  We choose a rectangular grid of complex numbers.  We iterate from each of the numbers until one of two things happen:
 • The iterate exceeds max(2,|c|) in absolute value
     - in which case we shade the seed according to how quickly it reached the limit.
 • We reach some pre-specified bailout number of iterates
     - in which case we assume the seed is in (Subscript[K, c]) and color it black

The same basic algorithm works for any polynomial, though there is a bit of work to figure out the appropriate bound.

Classification of fixed points and periodic orbits

The classification

Recall that a fixed point of a holomorphic function f is a point (Subscript[z, 0])∈ℂ such that f((Subscript[z, 0]))=(Subscript[z, 0]).  We classify fixed points as

Attracting if |f'((Subscript[z, 0]))|<1

Super-attracting if |f'((Subscript[z, 0]))|=0

Repelling if \[LeftBracketingBar]f'((Subscript[z, 0]))\[RightBracketingBar]>1

Neutral if |f'((Subscript[z, 0]))|=1

It’s worth mentioning that this same classification holds for the iteration of smooth real functions.

The theorems

The following theorems justify this language:

Theorem (The attracting fixed point theorem):
Suppose that (Subscript[z, 0]) is an attracting fixed point for the holomorphic function f.  Then, there is an r>0 such that (f^n)(z)->(Subscript[z, 0]) for all z∈(Subscript[D, r])((Subscript[z, 0]))={z:|z-(Subscript[z, 0])|<r}.

Proof:
Since |f'((Subscript[z, 0]))|<1, we may choose a positive λ<1 and an r>0 so that

|((f(z)-f((Subscript[z, 0])))/(z-(Subscript[z, 0])))|=((|f(z)-f((Subscript[z, 0]))|)/(|z-(Subscript[z, 0])|))<λ<1

for all z∈(Subscript[D, r])((Subscript[z, 0])).  It follows that

|f(z)-(Subscript[z, 0])|=|f(z)-f((Subscript[z, 0]))|<λ|z-(Subscript[z, 0])|

and, by induction, that

|(f^n)(z)-(Subscript[z, 0])|<(λ^n)|z-(Subscript[z, 0])|.

Since λ<1, we must have (f^n)(z)->(Subscript[z, 0]).□

Note that we have a bit more than just (f^n)(z)->(Subscript[z, 0]) - we know the rate at which it happens, namely, the convergence is similar to the convergence of ((f'((Subscript[z, 0])))^n).  In the case of a super-attracting orbit, we have the convergence is faster than (r^n) for every positive r.  There is a repelling fixed point theorem whose statement and proof is very similar.

Periodic points

Recall that (Subscript[z, 0]) is periodic if (f^n)((Subscript[z, 0]))=(Subscript[z, 0]) for some positive integer n.  If n is the least such integer, then the orbit of (Subscript[z, 0]) is just ((Subscript[z, 0]),(Subscript[z, 1]),(Subscript[z, 2]),…,(Subscript[z, n-1])) where f((Subscript[z, k]))=(Subscript[z, k+1]) for each k=0,1,…,n-2 and f((Subscript[z, n-1]))=(Subscript[z, 0]).

Suppose that (Subscript[z, 0]) has period n under iteration of f and let F=(f^n), so that (Subscript[z, 0]) is a fixed point of F.  We say that the orbit of (Subscript[z, 0]) is attractive orbit for f if (Subscript[z, 0]) is an attractive fixed point of F.  The notions of super-attracting, repelling, and neutral can be extended just as easily.

Computationally, this classification can be expressed entirely in terms of f using the chain rule.  For simplicity, suppose that (Subscript[z, 0]) has period 2 for f and define F=f∘f.   Then,

F'(z)=(d/(dz))(f(f(z))=f'(f(z))f'(z).

Thus, if f((Subscript[z, 0]))=(Subscript[z, 1]) and f((Subscript[z, 1]))=(Subscript[z, 0]), we have

F'((Subscript[z, 0]))=f'(f((Subscript[z, 0])))f'((Subscript[z, 0]))=f'((Subscript[z, 1]))f'((Subscript[z, 0])).

More generally, if (Subscript[z, 0]) has period n for f and F=(f^n), then

F'((Subscript[z, 0]))=f'((Subscript[z, 0]))f'((Subscript[z, 1]))(⋯f)'((Subscript[z, n-1])).

Thus, we can classify the orbit by multiplying the values of the derivatives of f at the points in the orbit; the function F need not be computed explicitly.

The Mandelbrot set

Visualizing the set

Here’s how to use the MandelbrotSetPlot command.

(* The simplest whole set *)
MandelbrotSetPlot[]

(* Setting the range and a few options *)
MandelbrotSetPlot[{-2-0.1I,-1.3+0.1I},
FrameTicks -> None, ImageSize->600,
ImageResolution->1200, MaxIterations->200,
ColorFunction->"StarryNightColors"]

(* A deeper zoom *)
MandelbrotSetPlot[{-1.52219977608+0.000107683187I,-1.5221966361+0.000104917046I},
MaxIterations->5000, FrameTicks -> False, ImageResolution->1200,
ColorFunction->"SunsetColors", ImageSize->600]

The critical orbit for (z^2)+c

A critical point of a function f:ℂ->ℂ is the same as in calculus 1 - i.e. a complex number (Subscript[z, 0]) such that f'((Subscript[z, 0]))=0.  As we’ll learn a bit later, complex dynamics are dominated by the behavior of the critical points.  The orbit of such a point is often called a critical orbit.

For the quadratic family (Subscript[f, c])(z)=(z^2)+c, we can understand the influence of the critical orbit fairly easily.  Of course, (Subscript[f, c]^′)(z)=2z so the origin is always the one critical point.  The independence of the critical point on the parameter is one of the appealing things about this family of function, by the way.

To understand the influence of the critical orbit, first recall the following theorem, which we proved above:

Theorem (The quadratic escape criterion):
Suppose that \[LeftBracketingBar]z\[RightBracketingBar]≥\[LeftBracketingBar]c\[RightBracketingBar] and \[LeftBracketingBar]z\[RightBracketingBar]>2. Then (Subscript[f, c]^n)(z)->∞ as n->∞.

Here’s a the key corollary:

Corollary:
If the critical orbit ever exceeds two in absolute value, then it diverges to ∞.

Proof: If \[LeftBracketingBar]c\[RightBracketingBar]>2, then \[LeftBracketingBar](Subscript[f, c])(0)\[RightBracketingBar]=\[LeftBracketingBar]c\[RightBracketingBar]>2.  Thus, the first iterate of the critical point is already larger than 2 in absolute value and the hypotheses of the escape criterion are satisfied so the critical orbit diverges.  If \[LeftBracketingBar]c\[RightBracketingBar]≤2 and any iterate ever exceeds 2 in absolute value then, again, the hypotheses are satisfied.\[Square]

Combined with our observations on totally disconnected Julia sets, this leads to the following refined and expanded version of the quadratic dichotomy:

If we iterate (Subscript[f, c]) from the critical point zero, then one of the following must happen:

The critical orbit stays bounded by 2 in absolute value, in which case the Julia set (Subscript[J, c]) is connected, or

The crucial orbit diverges to infinity, in which case the Julia set (Subscript[J, c]) is totally disconnected.

Thus, the orbit of the single critical point tells us a lot about the Julia set and the global dynamics of (Subscript[f, c]).

Definition

The Mandelbrot set, sometimes denoted just M, is the set of all complex numbers c such that the critical orbit of (Subscript[f, c])(z)=(z^2)+c stays bounded.  Equivalently, it is the set of all complex numbers c such that the Julia set of (Subscript[f, c]) is a connected set.

An algorithm

There is an algorithm to draw the Mandelbrot set that is very similar to the one for quadratic Julia sets.  We choose a rectangular grid of complex numbers representing c values - each of which determines a function (Subscript[f, c]).  For each c, we iterate (Subscript[f, c]) from the critical point zero until one of two things happen:
 • The iterate exceeds 2 in absolute value
     - in which case we shade the seed according to how quickly it reached the limit.
 • We reach some pre-specified bailout number of iterates
     - in which case we assume the seed is in M and color it black

Note the distinction - a Julia set is associated with a single function.  We iterate that single function from a collection of initial seeds and shade each seed according to the behavior.  The Mandelbrot set is a tool to help us understand an entire parametrized family of functions.  Each function is determined by a single complex number (the parameter).  For each choice of the parameter, we iterate the corresponding function from the critical point and shade the parameter according to the behavior.

The Julia set lives in the dynamical plane.

The Mandelbrot set lives in the parameter plane.

Parametrizing the components

As it turns out, we can get quite a solid handle on where these components come from.

The main cardioid

The largest and most conspicuous feature of the Mandelbrot set is the main cardioid.  It’s interior consists of all complex parameters c such that the function (Subscript[f, c])(z)=(z^2)+c has an attracting fixed point.  We can write down a pair of equations to describe the boundary of this domain:

(z^2)+c=z
2z=(e^(2π i t))

The first equation says we’re looking for a fixed point; the second says the fixed point should be neutral.

If we solve the second equation for z to get z=(e^(2π i t))/2, plug that into the first, and solve for c we get

c=(1/2)(e^(2π i t))-(1/4)(e^(4π i t))

Might looks a little tricky, but it’s easy to parametric plot it!

toPoint[z_]:={Re[z],Im[z]};
mainCardioid = ParametricPlot[
toPoint[Exp[2Pi*I*t]/2-Exp[4Pi*I*t]/4],
{t,0,1}]

The period two bulb

Clear[f,F,z,c];
f[z_] := (z^2)+c;
F[z_] = f[f[z]];
eqs = {F[z]==z,F'[z]==Exp[2Pi*I*t]};
Solve[eqs,{z,c}]

toPoint[z_]:={Re[z],Im[z]};
period2Bulb = ParametricPlot[
toPoint[-1+Exp[2Pi*I*t]/4],
{t,0,1}]

both = Show[{mainCardioid, period2Bulb}, PlotRange->All]

Show[{MandelbrotSetPlot[], both}]