Newton's Method (Complex)

We now show that Newton's method extends easily to one complex equation in one complex unknown. For example, let us see whether Newton can confirm that the roots of 

   z2-2z+2

are in fact 

   1+i and 1-i.

The Newton rule remains 

   z=z-f(z)/f'(z)

and so I hope it is obvious to you that a real starting guess will produce only real iterates and will not get us back to our complex roots. Please open my meandering diary and then run a few examples for yourself. Now, with a nonreal guess, our diary indeed has a happy ending. Now, if all we wanted were the roots, we would simply type 

   roots([1 -2 2])

at the matlab prompt. Our interest however, is not the destination but the journey. Along the way we will gain an appreciation for the limitation of Newton's method and we will grow our pallete of Matlab art tools.

In particular, we will dissect Isaac's attack on a polynomial by first assigning each root a color. Then, depending on where Isaac takes us we paint accordingly. For example, if Isaac takes us to the first root we paint it red, if he takes us to the other root we paint it blue, if Isaac loses his way we paint the point black. As he has a much better chance of losing his way in the complex plane we will produce a number of gorgeous Newton wastelands.

You will write code that dissects a user specfied quartic. To get started we dissect the particular 

    (z2-1)(z2+0.16)

Although the main Newton step is of course 

   z = z - (z^2-1)*(z^2+0.16)/(4*z^3 - 1.68*z)		(NS)

it is unclear which, if any, root this may lead to. It all depends on where we start. We shall see now that (NS) beautifully partitions the complex plane into 5 distinct regions - four so-called Newton Basins and one remaining Newton Wasteland. The Newton Basin associated with the root z=1 is defined to be all those starting points for which (NS) eventually produces z=1. The other basins are defined accordingly. The Newton Wasteland corresponds to all those starting points for which (NS) fails to converge to a root. For more details and examples please browse the Student Gallery and hit this Newton Basin site. With respect to our example, we choose the following color scheme 


   Basin of 1 is yellow
   Basin of -1 is red
   Basin of 0.4i is green
   Basin of -0.4i is blue
   Wasteland is black

With this scheme this pointilist driver produces the picture at right 

    for x = .15:.0025:.55,
    for y = -.15:.0025:.15,
        color = newt(x+i*y,1e-3,20);
        plot(x,y,color)
        hold on
    end
    end

Let us reverse engineer what newt must be up to.

Well, it takes a seed (x+i*y) and a tolerance (1e-3) and a maxiter (20) and returns a color ('y', 'r', 'b', 'g', or 'k') depending on where (NS) took the seed. The driver then paints that seed accordingly. The range of x and y specify the window of the complex plane to be painted, while the increment (0.0025) specifies the resolution. Here is the code that does the job.

Though this is indeed straightforward to understand, it does not exploit Matlab's ability to work directly on vectors and matrices. We shall see that the nested for loops can be contracted to a single (very fast) line. 

    x = .15:.0025:.55;
    y = -.15:.0025:.15;
    [X,Y] = meshgrid(x,y);
    Z = X+i*Y;

    for k=1:maxiter,

        Z = Z - (Z.^2-1).*(Z.^2+0.16)./(4*Z.^3 - 1.68*Z); 

    end

Here Z denotes the full grid of complex Newton seed's, and, thanks to the dot operator (recall quad) we may operate on all seed's simultaneously.

Now, upon completing this for loop, Z will be a matrix whose elements are either close to one of the four roots, or not. To find out which ones are close to 1 I recommend the find command 

    [i1,j1] = find(abs(Z-1)<0.1);

Now you simply place yellow at (x(j1),y(i1)). In order to see what is happening please consult this 9 seed diary.

If we paint the other roots accordingly, and leave the wasteland blank, we arrive at the finer (a markersize of 1) portrait below.

Our assignment this week is to produce such portraits for arbitrary quartics. Quartics are encoded by 5 complex scalars, for example, the quartic 

    2z4-4z3+(2-i)z2-iz+10

is simply 

    q = [2 -4 2-i -i 10]

In order to evaluate a polynomial at a given z one just types 

    f = polyval(q,z)

Now, the derivative of the polynomial above is simply 

    dq = [8 -12 4-2i -i]

And to evaluate it at a given z one just types 

    df = polyval(dq,z)

You may use polyval for this assignment but I ask that you write your own polyder. You should give it a new name and test it against polyder.

The last innovation for this assignment is the translation of the q into a string that will be used to title your portrait. Typing help strfun will lead to hours of string fun - after which you'll know all you need for this week. But just to be sure, the line 

    qlab = strcat(qlab,'+',num2str(q(k)),'z^',num2str(5-k))

ought to be helpful. Of course indiscriminate use of this will lead to titles like 

    1z4+0z3+-0.84z2+0z+-0.16z0

rather than the slick one that adorns our finer portrait. I know I can count on you to suppress leading 1s, kill terms with leading 0s, and never write +- or z0. You might use strrep to replace these unsightly strings.