Continuing my last blog entry, let us try the Newton Algorithm with two variables. First we define the function f(x,y). We use an example where we know the zeros exactly.

\(f(x,y) = (x+iy)^3-1.\)

The zeros are the three roots of unity. In EMT, we define the real and the imaginary parts of this function as f1(x,y) and f2(x,y).

>fe &= z^3-1 with z=x+I*y

                                      3
                             (I y + x)  - 1

>function f1(x,y) &= realpart(fe)

                                  2    3
                           - 3 x y  + x  - 1

>function f2(x,y) &= imagpart(fe)

                                 2      3
                              3 x  y - y

We can now plot the function.

>plot2d("abs(f1(x,y)+f2(x,y)*I)",r=1.2,>contour,>hue,>spectral,n=100);

We can now define the function f([x,y]) as a hybrid function in EMT which can be used as f(x,y) or f(v) numerically, and likewise the derivative matrix of this function, courtesy by Maxima. Furthermore, we define the Newton step.

\(f(v) = v – Df(v)^{-1} \cdot f(v).\)

>function f([x,y]) &= [f1(x,y),f2(x,y)]

                            2    3         2      3
                    [- 3 x y  + x  - 1, 3 x  y - y ]

>function Df([x,y]) &= jacobian(f(x,y),[x,y])

                      [    2      2              ]
                      [ 3 x  - 3 y     - 6 x y   ]
                      [                          ]
                      [                 2      2 ]
                      [    6 x y     3 x  - 3 y  ]

>function newtonstep ([x,y]) &= transpose([x,y]-invert(Df(x,y)).f(x,y));

The function is the following.

\begin{pmatrix}\frac{2\,x\,y^4+4\,x^3\,y^2-y^2+2\,x^5+x^2}{3\, \left(y^2+x^2\right)^2} & \frac{2\,y\,\left(y^4+2\,x^2\,y^2+x^4-x \right)}{3\,\left(y^2+x^2\right)^2} \\ \end{pmatrix}

To get the Latex code for this, use the tex() command of Maxima.

>tex(&factor(newtonstep(x,y)))
 \begin{pmatrix}\frac{2\,x\,y^4+4\,x^3\,y^2-y^2+2\,x^5+x^2}{3\,
  \left(y^2+x^2\right)^2} & \frac{2\,y\,\left(y^4+2\,x^2\,y^2+x^4-x
  \right)}{3\,\left(y^2+x^2\right)^2} \\ \end{pmatrix}

Let us test.

>v=[1.5,0.4]; loop 1 to 6; v=newtonstep(v), end;
 [1.11995,  0.197797]
 [0.988758,  0.043586]
 [0.998147,  -0.000885968]
 [1,  3.29448e-006]
 [1,  1.74668e-011]
 [1,  0]

Okay, [1,0] is one of the zeros. Finally, we plot one path of the Newton Algorithm.

>v=[0.5,0.3]; w=v; loop 1 to 10; v=newtonstep(v); w=w_v; end; w,
           0.5           0.3 
      0.794694     -0.665052 
      0.584502     -0.137812 
       1.21661      0.321028 
      0.994208      0.110139 
      0.987868   0.000509852 
       1.00015 -1.26776e-005 
             1 -3.78536e-009 
             1             0 
             1             0 
             1             0 
>plot2d(w'[1],w'[2],color=red,thickness=2,>add):

We see that the algorithm can take quite erratic paths. So, it is interesting from which points the algorithm converges to which zeros. We simulate that.

>function map res(x,y) ...
$v=[x,y];
$loop 1 to 100;
$   if all(v~=0) then return 1; endif;
$   vnew=newtonstep(v);
$   if all(vnew~=v) then
$       if vnew[1]>0 then return 3;
$       elseif vnew[2]>0 then return 4;
$       else return 5;
$       endif;
$   endif;
$   v=vnew;
$end;
$return 1;
$endfunction
>res(0.5,0.3)
 3

Since we used the „map“ keyword we can apply the code to many pairs [x,y]. It will vectorize to a matrix of X and Y values too. We use it to compute the results of 500×500 points. On my computer, this takes a few seconds.

>X=linspace(-1.2,1.2,500); Y=X';
>R=res(X,Y);
>plotrgb(R):

As you see there are three main regions of convergence. The image is a nice fractal between these main regions.