Computing the Julia Set and the Mandelbrot Set

The Julia Set

In this exercise, we will display some Julia sets and then compute the Mandelbrot set. Our ultimate goal is to understand why the algorithm we present for computing the Mandelbrot set works.

Recall  the definition of the filled Julia set. We are supposed to fix a complex number  c  and study the dynamics of the map  f(z) = (z^2) + c.  As we will soon see , all large values of  z  go to infinity. The filled Julia set is the collection of z's which do not go to infinity. There certainly exist such points, because the equation (f^n)(z) = z can always be solved over the complex numbers by the fundamental theorem of algebra, and therefore there exist n-cycles which do not go to infinity. 

Below is a Mathematica function. This  function defines an empty list ``mylist'' and then adds all points in a grid to the list. At the end, these points are plotted. Examine the code carefully until you understand why it works.

M1[incre_,size_]:=Module[{u,v,mylist},mylist={};For[u=-size,u≤size,u=u+incre,For[v=-size,v≤size,v=v+incre,mylist=Append[mylist,{u,v}];];];ListPlot[mylist,AspectRatio->Automatic]];

M1[.1, 2]

We are going to modify this program so it only draws grid points when they belong to the Julia set. Begin by defining our standard function. Please execute the following line.

f[z_] := z^2 + c;

Now a little lemma. Suppose both inequalities below are true: 

	| z |  >  | c |
	| z |  >  3
	
Then | f(z) | = | (z^2) + c | ≥  (((| z |))^(2 ))- | c |  ≥ (((| z |))^2) - | z | = | z | ( | z | - 1) ≥ | z | 2.  So for such points, f(z) is further from the origin than  z, and each additional iteration sends it yet further away. Consequently, (f^n)(z) -> ∞.

Today we will always deal with  c  satisfying  | c | < 3  So as soon as z is further than 3 away from the origin, it goes to infinity.

Consider then the following code. It cycles through all gridpoints exactly as the previous code did. But this time, each time the code comes to a grid point  u + i v, it sets  z = u + i v and then iterates to see if the point gets beyond 3 from the origin. The code only adds a gridpoint to the list if the point did not get beyond 3.  So this code should draw all points in the filled Julia set.

The new code we  add is set off from the original code above by extra blank lines.

M2[iterations_,incre_,size_]:=Module[{u,v,z,mylist},mylist={};For[u=-size,u≤size,u=u+incre,For[v=-size,v≤size,v=v+incre,z=u+I v;i=1;While[i<iterations&&Abs[z]≤3,z=f[z];i=i+1;];If[i==iterations,mylist=Append[mylist,{u,v}];];];];ListPlot[mylist,AspectRatio->Automatic]];

Try this with c = -1. Consider the line below. Notice that by changing  .1 to .02 we get a clearer picture. We can improve the accurate of this picture by increasing the maximum number of iterations from 10 to 100. But I recommend sticking with 10 and .1 if you don't want to waste a lot of time.

c = -1; M2[100, .1, 2]

c = -1; M2[100, .02, 2]

c = -1; M2[100, .01, 2]

c = .2 + .3 I; M2[100, .1, 2]

c = .2 + .3 I; M2[100, .02, 2]

Please compute the Julia sets corresponding to c = .5 + .2 i, c = -1.3,  c = -1.75, and c = i. 

In particular, thinking carefully about c =   .5 + .2i. The Julia set for this  c  is a Cantor set. To visualize this set, you may need to decrease 1 to .04. You'll still see some big connected blobs. You'll need to increase 10 to a larger value to notice that these blobs also break up.

When a calculation requires a very large number of simple arithmetic operations, it usually runs slowly in Mathematica and much faster in C, Java, Pascal, or Fortran. If you can, try the experiments above with one of the faster languages.

Later in these notes we'll study the Julia set of (z^2) + (.7 + .8 i). Here it is:

c =1.4 + 1.5 I; M2[5, .04, 2.1]

The Mandelbrot Set

The function  f(z) = (z^2) + c has derivative zero exactly when z = 0. For this reason, z = 0 is a very special point. We'll spend two or three days trying to understand the reason it is so special. Let me tell you in advance the three main theorems which we'll partly understand:

THEOREM 1: If f(z) = (z^2) + c  has an attractive cycle, then z = 0 is attracted to this cycle. (In particular, f(z) can have at most one attractive cycle.)

THEOREM 2: If z = 0 does not go to infinity, then the filled Julia set is connected. It may happen that this filled Julia set is one-dimensional with no interior, but these cases are rare and usually there is an open interior with interesting dynamics. which attracts z = 0.

THEOREM 2: If z = 0 goes to infinity, then the filled Julia set is just a Cantor set. Consequently there is no open interior set with interesting dynamics.


Using these theorems, we  give the rigorous definition of the Mandelbrot set. It is defined to be the set of all c such that under iteration by f(z) = (z^2) + c, the point z = 0 does not go to infinity. Consequently, it is the set of all c such that the filled Julia set is connected.

We can modify our inequalities above to prove that if  c  is a point of the Mandelbrot set, then | c | ≤ 2. We can also prove that if z = 0  goes to a point further than 2 from the origin, then it will go to infinity. You'll find this slight modification in the notes  handed out on Friday.

The code below draws our standard grid. Modify this code so it draws the Mandelbrot set. At the start of the inside piece, we have set c to the grid point. You should then start with z = 0 and iterate f(z) = (z^2) + c many times, testing to see if the result gets further than 2 from the origin. If it does, the point (u, v) should not be added to the list, but otherwise it belongs to the Mandelbrot set and should be added.

M[iterations_,incre_]:=Module[{u,v,c,mylist},mylist={};For[u=-2,u≤1,u=u+incre,For[v=-1.5`,v≤1.5`,v=v+incre,c=u+I v;i=1;If[0==0,mylist=Append[mylist,{u,v}];];];];ListPlot[mylist,AspectRatio->Automatic]];

M[20, .1]

Why The Julia Set Is A Cantor Set for Large c (Part 1)

During the rest of the worksheet, suppose that | c | > 2. We are going to explain why the Julia set looks like a Cantor set. As the argument proceeds, you should notice two things. First, notice that the argument is similar to the argument in class which gave us a Cantor set associated with (x^(2 ))+ c as soon as c < -2. Second, notice the special role played by z = 0 in the argument.

We first claim that when |z| > |c|, then z will go to infinity. Indeed, |f(z)| = |(z^2) + c| ≥|(z^2)| - |c| ≥ |(z^2)| - |z| = |z| (|z| - 1) > |z|,
so further iterates get larger and larger. Therefore, the points which remain bounded under iteration all belong to the disk of radius |c|. 

Now fix a c. Our first task is to draw this disk of radius |c|. Only points within this disk remain bounded. (Even some of these points go to infinity, as we'll see later.) I'll take c = (.7 + .8 i), but feel free to try some other c instead. Please execute the line below, or change it and then execute it.

c = 1.4 + 1.5 I; Abs[c]

We first want to draw the circle of radius |c| because we know that all points outside this circle will to go infinity. The equation of the circle is | z | = | c |. If we expand this condition  using z = x + i y and c = c1 + i c2, we get the equation

     (x^2) + (y^2) = (c1^2) + (c2^2)

To draw solutions of this equation, we could simply ask Mathematica to draw a circle.  But later we need to do the same thing for more complicated euqations. There is a better way. 

Mathematica can draw curves which solve f(x, y) = 0 where f is any reasonable function with two variables. Let's  try a simple example. We'll plot the ellipse   (x^2) + 2 (y^2) = 1. Below is the command to do it. Mathematica wants to draw the points where f(x, y) = 0, so we rewrite our equation with (x^2) + 2 (y^2) - 1 on one side and zero on the other side. The command ``Contours->{0}'' says to draw only the curve where f = 0, the command ``ContourShading->False'' says not to blacken the interior, and the command ``PlotPoints->50'' says to calculate more than the usual number of points for a more accurate drawing.

ContourPlot[x^2 + 2 y^2 - 1, {x, -2, 2}, {y, -2, 2},
	Contours->{0}, ContourShading->False,
PlotPoints->50]

Now we are ready to draw the points which truly interest us. We want to draw all z such that |z| = |c|. Below is the appropriate command.

Fig0 = ContourPlot[Abs[x + I y] - Abs[c], {x, -2.1, 2.1}, {y, -2.1, 2.1},
	Contours->{0}, ContourShading->False, PlotPoints->50]

Why The Julia Set Is A Cantor Set for Large c (Part 2)

And now the main point. Certainly points outside this circle go to infinity. But some points inside the circle may go outside in the next iterations. Those points will ultimately go to infinity. To find out which points inside leave after one iteration, we will find those points inside the circle which map exactly to the boundary of the circle under f. Notice (and this point is very important) that z = 0 is one of these points because (0^2) + c is on the boundary of the circle. We can find all the solutions by solving | f(z, c) | = |c|, which are do using the following command:

Fig1 = ContourPlot[Abs[f[x + I y]] - Abs[c], {x, -2.1, 2.1}, {y, -2.1, 2.1},
	Contours->{0}, ContourShading->False, PlotPoints->50]

The points exactly on this figure eight map exactly to the boundary of the disk of radius c after one iteration. The points outside this figure eight curve map outside the disk of radius  c  and consequently go to infinity. So only points inside the figure eight remain inside the disk after one iteration.

However, points inside the figure eight might still map outside the disk under two iterations. To see which points remain, we need to look for points which exactly map to the boundary of the disk under two iterations. These points are given by solving | (f^2)(z, c) | = |c| as below:

Fig2 = ContourPlot[Abs[f[f[x + I y]]] - Abs[c], {x, -2.1, 2.1}, {y, -2.1, 2.1},
	Contours->{0}, ContourShading->False, PlotPoints->50]

Only points inside these two figure eights will remain inside the disk after two iterations. However, some points inside the two figure eights may map to the boundary of the disk after three iterations. To discover them, we must compute

Fig3 = ContourPlot[Abs[Nest[f,x + I y,3]] - Abs[c], {x, -2.1, 2.1}, {y, -2.1, 2.1},
	Contours->{0}, ContourShading->False, PlotPoints->100]

Fig4 = ContourPlot[Abs[Nest[f,x + I y,4]] - Abs[c], {x, -2.1, 2.1}, {y, -2.1, 2.1},
	Contours->{0}, ContourShading->False, PlotPoints->100]

Let us put all of these diagrams together:

Show[Fig0, Fig1, Fig2, Fig3, Fig4]

To summarize, we are studying the map f(z) = (z^2) + (1.4 + 1.5 i). All points outside the disk of radius  | c |  go to infinity. However, some points inside this disk are mapped exactly to the boundary of the disk by f(z). The set of such points is a figure eight including the origin. Points outside this figure eight are mapped outside the disk and ultimately to infinity. So only points inside this figure eight remain bounded.

But  some points inside this figure eight map map to the boundary of the disk under two iterations. These points are exactly points on the two smaller figure eights. They map to the large figure eight under f, and so to the boundary of the disk under (f^2). Points outside these smaller figure eights ultimately map to infinity. The Julia set contains only points inside these figure eight curves.

Etc.

Please try this calculation with other values of c. Compare the resulting pictures to the pictures of Julia sets computed earlier, or by another more efficient Julia program.