# 1.12 Lesson 12

```# Lesson 12.  Newton's method in several variables
# =======================================
#
> restart: with(linalg):
Warning: new definition for   norm
Warning: new definition for   trace

# We have used Newton's method to solve equations in one
# variable.  Suppose now we wish to solve a system of m equations in m
# variables:
# f[1](x[1], ..., x[m])=0, ...  f[m](x[1],...,x[m]) = 0
#
# The basic idea is the same:
# - Consider linear approximations to the functions  near the initial
# guess.
# - Find a point where the linear approximations are all 0.
# - Use that as the next guess.  Repeat this process until it appears to
# converge.
#
> f1:= (x,y,z) -> x^2 + x*y*z + y^2*z^3 - 3;

2            2  3
f1 := (x,y,z) -> x  + x y z + y  z  - 3
> f2:= (x,y,z) -> x*y^2 - y*z^2 - 2*x^2;

2      2      2
f2 := (x,y,z) -> x y  - y z  - 2 x
> f3:= (x,y,z) -> x^2*y + y^2*z+ z^4-4;

2      2      4
f3 := (x,y,z) -> x  y + y  z + z  - 4
# We'd like to consider these as making up one mapping
# from R^3 to R^3. We could use Maple's "vector" structures
# to represent the points of R^3.
> V:= vector([a,b,c]);

V := [ a, b, c ]
# Note the absence of commas between the components.
# However,  I find it more convenient to use lists for this purpose.  If V
# = [x,y,z] is a list (representing a vector), op(V) is the sequence x,y,z,
# which can be used
# as the inputs to f1, f2 or f3.
> F:= V -> [f1(op(V)), f2(op(V)), f3(op(V))];

F := V -> [f1(op(V)), f2(op(V)), f3(op(V))]
> F([s,t,u]);

2            2  3         2      2      2   2      2      4
[s  + s t u + t  u  - 3, s t  - t u  - 2 s , s  t + t  u + u  - 4]
# We'll be using the list [x,y,z] a lot, so assign it to a name.
> xyz:= [x,y,z];

xyz := [x, y, z]
# Now the linear approximation to F(V) for V near V0 is
# F(V0) + J(V0) (V-V0) where J(V0) is the Jacobian matrix
# of F.  The "linalg" package has the "jacobian" command
# to produce it.  This takes two inputs, a list of expressions in the
# variables
# and a list of the variables.
> jacobian(F(xyz), xyz);

[                       3           2  2 ]
[ 2 x + y z  x z + 2 y z   x y + 3 y  z  ]
[                                        ]
[   2                  2                 ]
[  y  - 4 x   2 x y - z       - 2 y z    ]
[                                        ]
[              2              2      3   ]
[   2 x y     x  + 2 y z     y  + 4 z    ]
# We'd like to make this a function of [x,y,z].  Unfortunately Maple's
# matrices don't seem to work well with functions.  Most of the obvious
# approaches don't work.
> j:=unapply(",xyz);
Error, (in unapply)
variables must be unique and of type name

> j:= unapply(",(x,y,z));

j := (x,y,z) -> [tables and arrays not printed here]
> j(x,y,z);

[                       3           2  2 ]
[ 2 x + y z  x z + 2 y z   x y + 3 y  z  ]
[                                        ]
[   2                  2                 ]
[  y  - 4 x   2 x y - z       - 2 y z    ]
[                                        ]
[              2              2      3   ]
[   2 x y     x  + 2 y z     y  + 4 z    ]
> j(1,2,3);

[                       3           2  2 ]
[ 2 x + y z  x z + 2 y z   x y + 3 y  z  ]
[                                        ]
[   2                  2                 ]
[  y  - 4 x   2 x y - z       - 2 y z    ]
[                                        ]
[              2              2      3   ]
[   2 x y     x  + 2 y z     y  + 4 z    ]
# One approach that does work is to construct the function we want
# using lists of lists instead of matrices, and then use that with
# "matrix".  We take the expression for the Jacobian and convert it into
# a list of lists.
> convert(jacobian(F(xyz),xyz), listlist);

3           2  2     2                 2
[[2 x + y z, x z + 2 y z , x y + 3 y  z ], [y  - 4 x, 2 x y - z , - 2 y z],

2           2      3
[2 x y, x  + 2 y z, y  + 4 z ]]
# Now we make this list of lists into a function of x,y and z using
# "unapply".
> jlist:= unapply(",(x,y,z));

3           2  2
jlist := (x,y,z) -> [[2 x + y z, x z + 2 y z , x y + 3 y  z ],

2                 2                     2           2      3
[y  - 4 x, 2 x y - z , - 2 y z], [2 x y, x  + 2 y z, y  + 4 z ]]
# Now we can define a function that will take a list (representing a
# vector [x,y,z]) and return the Jacobian at that point as a matrix.
> j:= V -> matrix(jlist(op(V)));

j := V -> matrix(jlist(op(V)))
> j(xyz);

[                       3           2  2 ]
[ 2 x + y z  x z + 2 y z   x y + 3 y  z  ]
[                                        ]
[   2                  2                 ]
[  y  - 4 x   2 x y - z       - 2 y z    ]
[                                        ]
[              2              2      3   ]
[   2 x y     x  + 2 y z     y  + 4 z    ]
> j([1,2,3]);

[ 8  111  110 ]
[             ]
[ 0   -5  -12 ]
[             ]
[ 4   13  112 ]
# Now the linear approximation formula involves a multiplication
# (matrix)(vector).
# The usual "*" for multiplication of numbers (or expressions that
# represent numbers) can't be used for matrices and vectors, mainly
# because this is a noncommutative multiplication: A B and B A are
# different in general if A and B are matrices.  But Maple would assume
# they are the same.
> A*B - B*A;

0
# Instead, the multiplication symbol for use with matrices and vectors is
# "&*"
> A &* B - B &* A;

(A &* B) - (B &* A)
> matrix([[a,b],[c,d]]) &* [e,f];

[ a  b ]
[      ] &* [e, f]
[ c  d ]
# To actually perform linear algebra operations such as matrix
# multiplication (also addition, subtraction, and powers), you use the
# "evalm" command.
> evalm(");

[ a e + b f, c e + d f ]
# This produces a vector or matrix.  To make a vector into a list, we can
# use "convert".
> convert(",list);

[a e + b f, c e + d f]
# When I'm doing linear algebra, I find myself using the combination of
# evalm and convert(",list) so much that I define a single command to
# combine them:
> evl:= V -> convert(evalm(V),list);

evl := V -> convert(evalm(V), list)
# Here is the linear approximation to F(V) near V=[x0,y0,z0].
> evl(F([x0,y0,z0]) + j([x0,y0,z0]) &* (xyz-[x0,y0,z0]));

2                    2   3                                           3
[- x0  - 2 x0 y0 z0 - 4 y0  z0  - 3 + 2 x0 x + y0 z0 x + x0 z0 y + 2 y0 z0  y

2   2
+ x0 y0 z + 3 y0  z0  z,

2          2       2     2                            2
- 2 x0 y0  + 2 y0 z0  + 2 x0  + y0  x - 4 x0 x + 2 x0 y0 y - z0  y - 2 y0 z0 z,

2          2          4                     2                   2
- 2 x0  y0 - 2 y0  z0 - 3 z0  - 4 + 2 x0 y0 x + x0  y + 2 y0 z0 y + y0  z

3
+ 4 z0  z                                                           ]
> collect(",{x,y,z});

3                   2   2        2
[(2 x0 + y0 z0) x + (x0 z0 + 2 y0 z0 ) y + (x0 y0 + 3 y0  z0 ) z - x0

2   3
- 2 x0 y0 z0 - 4 y0  z0  - 3,

2                          2             2          2       2
(y0  - 4 x0) x + (2 x0 y0 - z0 ) y - 2 x0 y0  + 2 y0 z0  + 2 x0  - 2 y0 z0 z,

2                   2       3          2          2          4
2 x0 y0 x + (x0  + 2 y0 z0) y + (y0  + 4 z0 ) z - 2 x0  y0 - 2 y0  z0 - 3 z0

- 4                                                                     ]
# We'll make this a function of x0, y0 and z0.
> linapp:= unapply(",(x0,y0,z0));

linapp := (x0,y0,z0) -> [

3                   2   2        2
(2 x0 + y0 z0) x + (x0 z0 + 2 y0 z0 ) y + (x0 y0 + 3 y0  z0 ) z - x0
2   3
- 2 x0 y0 z0 - 4 y0  z0  - 3,

2                          2             2          2       2
(y0  - 4 x0) x + (2 x0 y0 - z0 ) y - 2 x0 y0  + 2 y0 z0  + 2 x0  - 2 y0 z0 z,

2                   2       3          2          2          4
2 x0 y0 x + (x0  + 2 y0 z0) y + (y0  + 4 z0 ) z - 2 x0  y0 - 2 y0  z0 - 3 z0
- 4                                                                     ]
> linapp(0,1,1);

[x + 2 y + 3 z - 7, x - y + 2 - 2 z, 2 y + 5 z - 9]
# Now we solve the equations from each component of this list to get
# the next Newton approximation.  For "solve" we need a set of
# expressions, not a list.
> solve({op(")},{x,y,z});

{z = 9/5, x = 8/5, y = 0}
# We'd like a command that will take the list [x0,y0,z0] as input and
# returns its result as a list [x,y,z].  We first use "op" on [x0,y0,z0] to
# get x0, y0, z0 which can be used as the input to "linapp", and we
# substitute the result of "solve" into [x,y,z].  One more thing: we want
# the result in floating point, not as an exact fraction (which after a few
# iterations could have very large numerator and denominator), so we
# use "evalf".
> newt:= V0 -> evalf(subs(solve({op(linapp(op(V0)))},{x,y,z}),[x,y,z]));

newt := V0 -> evalf(subs(solve({op(linapp(op(V0)))}, {x, y, z}), [x, y, z]))
> newt([0,1,1]);

[1.600000000, 0, 1.800000000]
# An alternative would be to use the "linsolve" command.
# linsolve(A,b) solves the vector equation A x = b, where A is a matrix
# and b is a vector.  In our case the equation to solve is F(V0) + j(V0) X
# = 0, i.e. j(V0) X = - F(V0), and we add this to V0 to get our next
# guess.  We can take the "-" outside the "linsolve".
> [0,1,1] - linsolve(j([0,1,1]),  F([0,1,1]));

[0, 1, 1] - [ -8/5, 1, -4/5 ]
# We want the vector addition evaluated, and the result returned as a
# list.
> newt2:= V0 -> evalf(evl(V0 - linsolve(j(V0),  F(V0))));\

newt2 := V0 -> evalf(evl(V0 - linsolve(j(V0), F(V0))))
> newt2([0,1,1]);

[1.600000000, 0, 1.800000000]
# According to my experiments, "newt" is considerably faster than
# "newt2".  One iteration of "newt2" more than 4 times as long as one
# iteration of "newt".  I'm not sure exactly why.
> V:= [0,1,1]: timer:= time():\
for jj from 1 to 20 do V:= newt(V) od: \
time()-timer;

.417
> V:= [0,1,1]: timer:= time():\
for jj from 1 to 20 do V:= newt2(V) od:\
time()-timer;

1.766
> V:= [0,1,1];

V := [0, 1, 1]
> for jj from 1 to 10 do V:= newt(V) od;

V := [1.600000000, 0, 1.800000000]

V := [-.4053571430, 2.380952381, 1.260183553]

V := [.1238168236, 1.962724782, 1.080486675]

V := [.3893816263, 2.074499230, .8377064543]

V := [.426964196, 1.979477834, .818450908]

V := [.4302423924, 1.978996584, .8151083308]

V := [.430244248, 1.978999685, .8150943555]

V := [.430244241, 1.978999694, .815094354]

V := [.430244244, 1.978999690, .8150943547]

V := [.430244243, 1.978999694, .8150943547]
# It's settled down nicely to what appears to be a solution.  The last
# digit or two changes at each iteration because of roundoff error.

> F(V);

-8         -8        -7
[.7*10  , -.12*10  , .10*10  ]
# Try a different starting point.
> V:= [1,0,0];

V := [1, 0, 0]
> V:= newt(V);

V := [x, y, z]
# The first step didn't seem to work.  What was the linear
# approximation?
> linapp(1,0,0);

[2 x - 4, - 4 x + 2, y - 4]
# Well, this clearly can't be solved for x, y and z.  We were unlucky to
# choose an initial V where the Jacobian is singular.  The "solve"
# command then produced nothing.  Curiously, if "subs" is given only
# one input it returns that input unchanged.
> subs([x,y,z]);

[x, y, z]
> subs(solve({op(linapp(1,0,0))},{x,y,z}),[x,y,z]);

[x, y, z]
# It would be better if "newt" detected this fact and produced an error
# message instead of returning [x,y,z].  We can do this by replacing
# "subs".  The replacement doesn't do anything differently, it just
# produces an error if there is only one input instead of two.
> newsubs:= proc(x,y)   subs(x,y) end;

newsubs := proc(x,y) subs(x,y) end

> newt:= subs(subs=newsubs, eval(newt));

newt := V0 -> evalf(newsubs(solve({op(linapp(op(V0)))}, {x, y, z}), [x, y, z]))
> newt([1,0,0]);
Error, (in newsubs) newsubs
uses a 2nd argument, y, which is missing

# Another starting point.
> V:= [0.2,1.3,1];

V := [.2, 1.3, 1]
> for jj from 1 to 10 do V:= newt(V) od;

V := [-1.071128893, 5.499445912, -.611168750]

V := [-1.049215594, 2.534205466, -.6192708733]

V := [-1.204204108, .5177684731, -.5997565973]

V := [-1.039618373, -7.863509003, -18.14022514]

V := [-14.91802573, -6.837306518, -13.63820635]

V := [-10.43664265, -5.919103138, -10.29934921]

V := [-7.234861711, -5.100300579, -7.801733768]

V := [-5.029357694, -4.362483323, -5.938344249]

V := [-3.537389653, -3.688691206, -4.553543125]

V := [-2.538975612, -3.064095347, -3.529896885]
# There's no sign yet of this one settling down.  Eventually it does, but
# to a different solution than the last one.
> for jj from 1 to 20 do V:= newt(V) od:
> V, newt(V);

[1.174658303, 1.674928536, .5655524760], [1.174658306, 1.674928535, .5655524741]
# How can we find out how many solutions there are?
#
# Our equations in this case involve polynomials, so "solve" might be
# useful.  Unfortunately, it isn't.
> # s:= solve({f1(x,y,z), f2(x,y,z), f3(x,y,z)}, {x,y,z});
# I put "#" before this one in the worksheet because I don't want
# you to try it: it seems to take forever.  Presumably the solution is
# extremely complicated.  But we can be a bit less ambitious and
# get some results.
# The intersections of f1(x,y,z)=0 and f2(x,y,z)=0 will correspond to
# the following curve in the yz plane.
> r1:= resultant(f1(x,y,z), f2(x,y,z), x);

2  4      3  5         2      4  6       2  3         4  3      5  4
r1 := y  z  - 4 y  z  + 12 y z  + 4 y  z  - 24 y  z  + 36 + y  z  + 2 y  z

3      6  3      4      3  4
- 6 y  z + y  z  - 3 y  + 2 y  z
# And similarly for f2 and f3.
> r2:=resultant(f1(x,y,z), f3(x,y,z),x);

5  3    3  6      3  2      3  3           2  5        4    4  2    6  6
r2 := y  z  + y  z  - 4 y  z  + 8 y  z  + 16 + 2 y  z  + 6 y z  + y  z  + y  z

7  3    8      4      3        5  4      4  3      2             2
- 2 z  y  + z  - 8 z  + 6 y  z - 2 y  z  - 6 y  z  + 9 y  - 24 y - 8 y  z
# We'll see the results next time.
>

```