#
#
# 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 ]
--------------------------------------------------------------------------------
