2.7 Newton's method for equations in many variables

```#
#
#                                              WORKSHEET#18
#
#                          Newton's method for equations in many variables
#
#
#
--------------------------------------------------------------------------------
# In this worksheet we use Newton's method to solve simultaneously equations of the form
#
#                                            f(x,y)=0 and g(x,y)=0.
#
# First we explore the linear algebra package, learning the syntax of some of the  vector and
# matrix commands.  In particular the product of matrices A, B is M=A&*B, or if we want
# to evaluate the product we use M=evalm(A&*B).
--------------------------------------------------------------------------------
> with(linalg):
> ?vector
--------------------------------------------------------------------------------
> v:=vector([3,4]);

v := [ 3, 4 ]
--------------------------------------------------------------------------------
> ?matrix
--------------------------------------------------------------------------------
> M:=matrix(2,2,[[1,2],[3,4]]);

[ 1  2 ]
M := [      ]
[ 3  4 ]
--------------------------------------------------------------------------------
> ?evalm
--------------------------------------------------------------------------------
> evalm(M&*v);

[ 11, 25 ]
--------------------------------------------------------------------------------
> evalm(v&*M);
Error, (in linalg[multiply]) matrix dimensions incompatible

--------------------------------------------------------------------------------
# Even though v is printed as a row vector, Maple is considering it as a column vector.  This
# is strange.
#
# Now consider the following 2 functions f(x,y) and g(x,y).  We want to find all (x,y)
# satisfying both f(x,y)=0 and g(x,y)=0.  We use the plots package to get an idea where
# there might be solutions and then the fsolve command to find them to better accuracy.
--------------------------------------------------------------------------------
> f:=exp(-x^2-y^2)-x^2;

2    2     2
f := exp(- x  - y ) - x
--------------------------------------------------------------------------------
> g:=cos(x*y)-y;

g := cos(x y) - y
--------------------------------------------------------------------------------
> with(plots):
--------------------------------------------------------------------------------
> graphf:=implicitplot(f,x=-2..2,y=-2..2):
--------------------------------------------------------------------------------
> graphg:=implicitplot(g,x=-4..4,y=-2..2):
--------------------------------------------------------------------------------
> display(graphf,title=`f(x,y)=0`);
--------------------------------------------------------------------------------
#

** Maple V Graphics **

--------------------------------------------------------------------------------
> display(graphg,title=`g(x,y)=0`);
--------------------------------------------------------------------------------
#

** Maple V Graphics **

--------------------------------------------------------------------------------
> display([graphf,graphg],title=`f(x,y)=0 and g(x,y)=0`);
--------------------------------------------------------------------------------
#

** Maple V Graphics **

--------------------------------------------------------------------------------
# There are 2 solutions indicated in the plot.  We use the mouse to point at the intersections
# and estimate their coordinates.  Notice also that the solutions seem to be symmetric about
# the y-axis.  This is evident from the equations themselves.
--------------------------------------------------------------------------------
> fsolve({f,g},{x,y},{x=-0.6..-0.5,y=0.8..0.9});

{x = -.5772352657, y = .8751057995}
--------------------------------------------------------------------------------
> fsolve({f,g},{x,y},{x=0.5..0.6,y=0.8..0.9});

{x = .5772352657, y = .8751057995}
--------------------------------------------------------------------------------
# To use Newton's method we have to make the two functions f(x,y), g(x,y)  into a single
# function F(x,y).  This is done by setting X equal to the 2-dimensional vector [x,y] and
# then defining F[x,y]=[f(x,y), g(x,y)].  The simultaneous equations f(x,y)=0 and  g(x,y)=0
# are then equivalent to the one equation F(x,y)=0.  Newton's method applied to this
# equation is the iteration U[k+1]=U[k]-K(U[k])&*F(U[k]), where U[0] is some initial
# approximation to the solution of F(x,y)=0 and K(U[k])  is the inverse of the jacobian
# matrix J evaluated at U[k].
--------------------------------------------------------------------------------
> X:=[x,y];F:=vector([f,g]);

X := [x, y]

2    2     2
F := [ exp(- x  - y ) - x , cos(x y) - y ]
--------------------------------------------------------------------------------
> ?jacobian
--------------------------------------------------------------------------------
> J:=jacobian(F,X);

[              2    2                      2    2  ]
[ - 2 x exp(- x  - y ) - 2 x  - 2 y exp(- x  - y ) ]
J := [                                                  ]
[        - sin(x y) y           - sin(x y) x - 1   ]
--------------------------------------------------------------------------------
> K:=inverse(J);

[                                    2    2      ]
[       sin(x y) x + 1      y exp(- x  - y )     ]
[ - 1/2 --------------      ----------------     ]
[             %1                   %1            ]
K := [                                                ]
[                                   2    2       ]
[        sin(x y) y       x (exp(- x  - y ) + 1) ]
[    1/2 ----------     - ---------------------- ]
[            %1                     %1           ]

2        2    2                      2    2              2
%1 := x  exp(- x  - y ) sin(x y) + x exp(- x  - y ) + sin(x y) x  + x

2        2    2
- y  exp(- x  - y ) sin(x y)
--------------------------------------------------------------------------------
> det(J);

2        2    2                        2    2                2
2 x  exp(- x  - y ) sin(x y) + 2 x exp(- x  - y ) + 2 sin(x y) x  + 2 x

2        2    2
- 2 y  exp(- x  - y ) sin(x y)
--------------------------------------------------------------------------------
# Notice that det(J)=2*(%1).  Substituting values for x,y into matrices whose entries are
# functions in x and y does seem to work directly.  We get around this difficulty by defining
# a substitution operator subU:w->subs(x=U[1],y=U[2],w), where w is an expression
# involving x and y, and then using the map command.
--------------------------------------------------------------------------------
> U:=vector([3,4]);

U := [ 3, 4 ]
--------------------------------------------------------------------------------
> sub(U,J);\

sub(U, J)
--------------------------------------------------------------------------------
> subU:=w->subs(x=U[1],y=U[2],w);

subU := w -> subs(x = U[1], y = U[2], w)
--------------------------------------------------------------------------------
> ?map
--------------------------------------------------------------------------------
#
> map(subU,J);

[ - 6 exp(-25) - 6    - 8 exp(-25)  ]
[                                   ]
[    - 4 sin(12)    - 3 sin(12) - 1 ]
--------------------------------------------------------------------------------
> evalf(");

[                              -9 ]
[ -6.000000000  -.1111035509*10   ]
[                                 ]
[  2.146291672      .609718754    ]
--------------------------------------------------------------------------------
> map(subU,F);

[ exp(-25) - 9, cos(12) - 4 ]
--------------------------------------------------------------------------------
> U:=vector([1.0,1.0]);

U := [ 1.0, 1.0 ]
--------------------------------------------------------------------------------
# The next command is just the loop for Newton's method.  Try executing the loop using
# the starting approximation U:=[1,1] instead of [1.0,1.0].
--------------------------------------------------------------------------------
> for k from 1 to 10 do U:=evalm(U-map(subU,K)&*map(subU,F)) od;

U := [ .6287374296, .9200144805 ]

U := [ .5768905349, .8778760621 ]

U := [ .5772347283, .8751061907 ]

U := [ .5772352656, .8751057995 ]

U := [ .5772352657, .8751057995 ]

U := [ .5772352657, .8751057995 ]

U := [ .5772352657, .8751057995 ]

U := [ .5772352657, .8751057995 ]

U := [ .5772352657, .8751057995 ]

U := [ .5772352657, .8751057995 ]
--------------------------------------------------------------------------------
> U:=vector([5.0,7.0]);

U := [ 5.0, 7.0 ]
--------------------------------------------------------------------------------
> for k from 1 to 15 do U:=evalm(U-map(subU,K)&*map(subU,F)) od;

U := [ 2.500000000, 20.49523077 ]

U := [ 1.249999999, 20.89909885 ]

U := [ .6249999993, 16.29303250 ]

U := [ .3124999996, -19.69774367 ]

U := [ .1562499998, -.17646234 ]

U := [ 2.314343954, .994163725 ]

U := [ 1.159802256, .6984130781 ]

U := [ .7094316789, .8173749252 ]

U := [ .5912471640, .8692845410 ]

U := [ .5774296414, .8750366550 ]

U := [ .5772353027, .8751057879 ]

U := [ .5772352656, .8751057995 ]

U := [ .5772352657, .8751057995 ]

U := [ .5772352657, .8751057995 ]

U := [ .5772352657, .8751057995 ]
--------------------------------------------------------------------------------
> U:=[6.0,-7.0];

U := [6.0, -7.0]
--------------------------------------------------------------------------------
> for k from 1 to 15 do U:=evalm(U-map(subU,K)&*map(subU,F)) od;

U := [ 2.999999999, -8.945943412 ]

U := [ 1.499999999, -20.15236713 ]

U := [ .7499999994, -17.42991070 ]

U := [ .3749999998, 16.30880845 ]

U := [ .1874999998, -.57847038 ]

U := [ 3.017101614, .845843852 ]

U := [ 1.508647469, .4844854263 ]

U := [ .8301097561, .7234015802 ]

U := [ .6229780443, .8502296093 ]

U := [ .5793360811, .8741805493 ]

U := [ .5772398149, .8751041367 ]

U := [ .5772352657, .8751057995 ]

U := [ .5772352657, .8751057995 ]

U := [ .5772352657, .8751057995 ]

U := [ .5772352657, .8751057995 ]
--------------------------------------------------------------------------------
> ?inverse
#
--------------------------------------------------------------------------------
> U:=vector([10.0,10.0]);

U := [ 10.0, 10.0 ]
--------------------------------------------------------------------------------
# In the next loop we have changed the order of substitution and inversion.  Previously we
# inverted the jacobian matrix symbolically (that is before substitution of values, when the
# entries were still expressions in x,y) and in the next loop we first substitute values into the
# jacobian J and then take the inverse numerically.  This is probably better.
--------------------------------------------------------------------------------
> for k from 1 to 15 do U:=evalm(U-inverse(map(subU,J))&*map(subU,F)) od;

U := [ 4.999999999, 18.47905426 ]

U := [ 2.500000000, 35.08047587 ]

U := [ 1.250000000, -95.35587003 ]

U := [ .6250000002, -26.02063569 ]

U := [ .3125000000, -10.30840157 ]

U := [ .1562500000, -1.348344440 ]

U := [ 2.543349544, .360638710 ]

U := [ 1.273633028, .5629845190 ]

U := [ .7546782886, .7713538881 ]

U := [ .6026015271, .8621271627 ]

U := [ .5778957510, .8748316379 ]

U := [ .5772357094, .8751056442 ]

U := [ .5772352656, .8751057995 ]

U := [ .5772352657, .8751057995 ]

U := [ .5772352657, .8751057995 ]
--------------------------------------------------------------------------------

```