Fractals and Newton's Method

In[32]:=

  Clear[f]
  f[z_] := z^4 - 1.0

In the complex plane, this function has four roots, z1 = 1, z2 = -1, z3 = i, and z4 = -i. Our error analysis guarantees that if we start near any of these roots, we will converge to the root. To see what happens for initial guesses far from any of the roots, we will use a color map. We associate a different color with each root and then color each initial guess in the plane according to the root to which it converges.

In[33]:=

  roots = {1, -1, I, -I}

Out[33]=

  {1, -1, I, -I}

We assign the values 0, 1, 2, and 3 to the four roots of the equation. Anything that is not a root will be assigned the value 4.

In[34]:=

  Clear[rootPosition]
  rootPosition[x_] := 0 /; Abs[x - roots[[1]]] < 10^-3
  rootPosition[x_] := 1 /; Abs[x - roots[[2]]] < 10^-3
  rootPosition[x_] := 2 /; Abs[x - roots[[3]]] < 10^-3
  rootPosition[x_] := 3 /; Abs[x - roots[[4]]] < 10^-3
  rootPosition[x_] := 4

We assign four colors to these four roots, using the RGB (red/green/blue) color system.

In[35]:=

  Clear[colorMap]
  colorMap[f_] := Which[ f < .2, RGBColor[.25,1,.25],
                         f < .4, RGBColor[1,0,1],
                         f < .6, RGBColor[0,1,1],
                         f < .8, RGBColor[0,0,1],
                         True,   RGBColor[1,1,1]]

Using a grid of initial values in [xMin,xMax] x [yMin,yMax], apply Newton's method to each initial guess and record the roots in a two dimensional table. We don't want our poor computer to have to find the derivative of f[z] for every initial guess, so we will do this once ahead of time. Applying Newton's method to f[z] is the same as finding a fixed point of the function,

In[36]:=

  Clear[g]
  g[z_] = N[Simplify[z - f[z]/f'[z]]]

Out[36]=

  0.25
  ---- + 0.75 z
    3
   z

The following command required 38 minutes on a Macintosh Quadra 700. In contrast, FindRoot required 22 hours for a 150 x 150 grid! FindRoot symbolically calculated the derivative of f 22,500 times. You may want to avoid evaluating the rest of the commands in this notebook. You can just scroll through the rest of the notebook to see the results.

In[37]:=

  Clear[x,y,z]
  xGridPoints = 200; yGridPoints = 200;
  xMin = -2.002; xMax = 1.998;
  yMin = -2.001; yMax = 1.999;
  deltaX = N[(xMax-xMin)/xGridPoints];
  deltaY = N[(yMax-yMin)/yGridPoints];
  table = Table[FixedPoint[g, x + y I],
  	{y,yMin,yMax,deltaY}, {x,xMin,xMax,deltaX}];//Timing

Out[37]=

  {2305.13 Second, Null}

Now apply rootPosition to the table. The result is saved to a file, so that if the system crashes later, our time-consuming calculations will not have been lost.

In[38]:=

  tableMap = Map[rootPosition, table,{2}];
  
  (* There must be a 0 and a 4 in tableMap in order *)
  (* for its values to scale to the correct colors. *)
  
  If[Max[tableMap] != 4,  tableMap[[1,1]] = 4];
  If[Min[tableMap] != 0,  tableMap[[1,2]] = 0];
  		
  (* Save the tableMap to a file named "tableMapnnn", *)
  (* where nnn is the value of xGridPoints.			*)
  
  Save[StringJoin["tableMap",ToString[xGridPoints]],tableMap]

Now let's show which initial points converged to which roots. (Open the next group of cells to see the graph.)

In[39]:=

  density = ListDensityPlot[tableMap, Mesh -> False,
    FrameTicks -> 
     {{{0,"-2"},{xGridPoints/4,"-1"},{xGridPoints/2,"0"},
  	{3 xGridPoints/4,"1"},{xGridPoints,"2"}},
  	{{0,"-2"},{yGridPoints/4,"-1"},{yGridPoints/2,"0"},
  	{3 yGridPoints/4,"1"},{yGridPoints,"2"}}},
    ColorFunction -> colorMap];

Out[39]

This is a fractal. If we zoom in on one of the lobes, at [0.6,1] x [0.6,1], the complexity of the image does not diminish. Let's take a look. We can define a function, newtonMap, to automate all of the necessary steps. The definition of newtonMap is contained in an initialization cell in the cell group containing the notebook title, at the beginning of this notebook.
The zoom region is offset slightly from [0.6,1] x [0.6,1] because, for the function f(z) = z^4 - 1, Newton's method will not converge with initial guesses that are real multiples of 1 + i.

Exercise 4

In[40]:=

  newtonMap[f, 0.602, 0.998, 0.601, 0.999, 200, 200];

Out[40]

We can't stop now! Let's zoom in one more time, on the region
[0.732, 0.744] x [0.768, 0.780]. Open the next group of cells to see the image.

In[41]:=

  newtonMap[f, 0.732, 0.744, 0.768, 0.780, 200, 200];

Out[41]

Exercise 5

Exercise 6