Math 171 - Alexiades
                          Lab 13
            Fractals from Newton iterations

---

It is interesting to examine the convergence of Newton iterations to roots for various starting points.
Theory tells us that for each root there is a region, "basin of attraction", such that if the initial guess
lies in this region then the Newton iteration will converge to the target root eventually.
When the function F(z) has several roots, each root has its own basin of attraction.
What happens near the borders of such regions belonging to different roots? Is there a clear demarkation between regions?
It is remarkable that not only there is no clear demarkation but in fact the "boundaries" of basins of attraction are fractal structures!
The goal of this lab is to produce such fractals out of Newton iterations (for the cubic F(z) = z³−1) .

Consider the cubic F(z) = z³ − 1. Its roots are the three
      cubic roots of unity:   root_m = e^{i (2π m/3)},   m = 0, 1, 2   i.e.   1 , e^{i (2π/3)} , e^{i (4π/3)}
which, by Euler's formula (e^iθ=cosθ+i sinθ), are:   1 , cos(2π/3) ± i sin(2π/3) = ( −1 ± i√3 ) / 2 .

A. To do:

Set the three roots:   root1 = 1+0i;   root2 = exp(i*2*pi*1/3);   root3 = exp(i*2*pi*2/3);

Pick maxIT = 20 , radius of convergence R = 1 , (TOL not needed, accuracy is of no concern here).

Consider a window: [−1 , 1] × [−1 , 1] on the xy-plane.

For integer M , construct a grid of (M+1) × (M+1) z-grid points:   z_ij = x_i + i y_j  (as complex numbers).

Initialize an array for colors:   Colors = zeros(M+1, M+1)

Take each z-grid point as starting point (z₀) and compute its orbit under the Newton iteration:
      z_n+1 = z_n − F(z_n) / F'(z_n) ,   n = 0:maxIT   ( after calling [Fn, DFn] = FCN( zn ) , of course ).

For each iterate z_n:  if z_n is within radius R of root_m , m=1,2,3 then set color for the starting z_ij to be m:   Colors(jy,ix) = m;

Once all the z-grid points have been processed and assigned a color,
  plot the Colors array:   pcolor( Colors ), shading('flat') .   Note: No need to rescale the Colors array, only using 3 colors here.

B. Prepare for submital:
  1. Write a Matlab function NewtFrac( M ) that implements the above.
  2. With M=100, fractal structure will be visible (good for debugging).
    With M=200, fractal will look nicer, but it will take longer to run...
    With M=400, fractal will look much nicer, but it will take much longer to run...
  3. Save the plot with M=200 into a file named "NewtFrac_200.jpg"
  4. Play with maxIT and R to see the effect (if any).
  5. Clean up your NewtFrac.m code to contain ONLY what's needed (in particular, no interactive input).
  6. Publish it:   publish( NewtFrac(200) , 'pdf' )
    If you clean up your code well enough, it will work. View the PDF to make sure it contains both the code and the plot.

C. Submit on Canvas:
  Upload just the NewtFrac.pdf
  If you can't make 'publish' work, upload your "NewtFrac.m" code (as plain text file), and also upload your "NewtFrac_200.jpg" file.