% times noreferences

\begin{center}
\bf\Large Solving systems of nonlinear equations
\end{center}

\vspace{0.2in}

This handout describes Newton's method for solving a system of nonlinear
equations.  In essence it works by converting the problem into the problem
of solving some systems of linear equations.

There is an accompanying lab project, which is to be done in groups of 
three, and which is due March 15 (right before spring break).  I will pass 
around a sign-up sheet for the groups.

Before reading further, review the handouts and your class notes for the
material on linear and affine mappings.  Also note that the various pieces
of Maple code in this document are available electronically.  You can
access them by going to the FileViewer, choosing MathLabNotes, then
Math314S2Notes, and then newt.ms.  If you don't know how to do this, 
ask the lab attendant.  You may copy electronically from the given worksheet
to your worksheet.  In fact this is far better than copying by hand, which
would in all likelihood introduce errors.  

\midhead{Approximation by affine mappings} 

Suppose we have a mapping \mp[[ f || \R^n || \R^k ]] which is not affine.  Our 
objective is to approximate $f$ by affine mappings.  This is very much the 
point of calculus!

Let's start with a differentiable mapping \mp[[ f || \R || \R ]].  Fix 
$x_0 \in \R$.  The {\it best affine approximation to $f$ at $x_0$\/} is the 
affine mapping \mp[[ Lf || \R || \R ]] given by 
$$Lf(x) = f(x_0) + f'(x_0) (x - x_0).$$%
(The notation ``$Lf$'' is just the name of this mapping.)
Note that $Lf$ is affine.  (Are you sure you understand why?)  In fact,
$Lf(x)$ is characterized by the following properties:
\begin{itemize}
\item $Lf$ is affine
\item $Lf(x_0) = f(x_0)$
\item $(Lf)'(x_0) = f'(x_0)$.
\end{itemize}
Otherwise said, $Lf$ is the unique affine mapping which agrees up to first 
order with $f$ (at $x_0$).  

For example, let $f(x) = \sin(x)$, and let $x_0 = 1.3$.  The graphs of 
$y = f(x)$ and $y = Lf(x)$ are exhibited below:

\widepost{160}{80}{
{150 NorthSouthSegment newpath -50 0 moveto 150 0 lineto stroke
/PI {3.14159} def /x0 {1.3} def /todegrees {180 mul PI div} def
{newpath 0 0 moveto 0 x0 todegrees sin 20 mul lineto AltDotStroke
{(1.3) t-r} LabelBelow} x0 20 mul xput
{20 div todegrees sin 20 mul} -16 142 1 FunctionArrowDraw
{20 div x0 sub x0 todegrees cos mul
x0 todegrees sin add 20 mul} -16 142 1 FunctionArrowDraw
{{{(y = f(x)) t-r} LabelRight} 142 20 div todegrees sin 20 mul yput
{{(y = Lf(x)) t-r} LabelRight} 142 20 div x0 sub x0 todegrees cos mul
x0 todegrees sin add 20 mul yput} 142 xput} -50 xput}

Notice that the graph of $y = Lf(x)$ is just the tangent line to the 
graph of $y = f(x)$ at the point where $x = 1.3$.

Now suppose we have a mapping \mp[[ f || \R^2 || \R ]].  For such mappings,
we will always assume that both partial derivatives exist.  Fix 
$(x_0, y_0) \in \R^2$.  Then the {\it best affine approximation to $f$ at 
$(x_0,y_0)$\/} is the affine mapping \mp[[ Lf || \R^2 || \R ]] given by 
$$Lf(x,y)\ =\ f(x_0,y_0) + {\pa f \over \pa x}(x_0, y_0)(x-x_0)
  + {\pa f \over \pa y}(x_0, y_0)(y - y_0).$$%
Notice that (e.g.) we are computing the partial derivative of $f$ 
\WRT\ $x$ at the point $(x_0,y_0)$, and then multiplying that by $x-x_0$.
For example, suppose that $f(x,y) = x^3 + 3 x^2 y + y^3 + 2$.  Then 
$$Lf(x, y)\ =\ (x_0^3 + 3 x_0^2 y_0 + y_0^3 + 2) +
  (3x_0^2 + 6x_0y_0)(x-x_0) + (3x_0^2 + 3y_0^2)(y-y_0).$$

\midhead{On solving equations}

Given a system of linear equations, we can solve for and describe the 
solution set exactly, using Gauss-Jordan elimination.

For a system of nonlinear equations, there is in general no way to describe
the solution set exactly.  Even when one has a single equation, involving
only one variable, one cannot in general get an exact description of the
solution set.

Even describing the solution set approximately is going to be a horrendous
problem (in general) if there are infinitely many solutions.  So we'll
restrict attention to systems of nonlinear equations which we think 
{\it probably\/} have only finitely many solutions.  More specifically,
assume that {\it the number of variables equals the number of equations}.
(What additional condition is needed in the {\it linear\/} case to insure 
that there are only finitely many solutions?)

In the next section we describe Newton's method for finding approximate
solutions to systems of nonlinear equations.  To use the method, you
first have to come up with an initial guess.  This need not be close to
an actual solution, and may more or less be chosen at random.  But whether
or not the method works, and which solution it ultimately finds will depend
on your initial guess.  So you have to try many initial guesses.  But it
is not easy to know how many initial guesses should be tried, since there
is no easy way to tell how many solutions the system has.  The whole
process is much like hunting for Easter eggs: you never know quite where to
look, and you never know for sure that you've found them all.

\midhead{Newton's method}

Now we illustrate how the idea of best affine approximation is useful for
finding approximate solutions to systems of nonlinear equations.  Consider
the system:
\begin{eqnarray*}
x^3 + 3 x^2 y + y^3 + 2 & = & 0\\
\sin(x) + \sin(y) + x y + 5 & = & 0.
\end{eqnarray*}
Our objective is to find points $(x,y) \in \R^2$ which are approximate
solutions to both equations.  The method is exactly Newton's method,
generalized from the one variable case you have seen in calculus.  It goes like 
this.  First pick a point $(x_0, y_0) \in \R^2$, more or less at random.  
(This point is called the {\it initial guess}.)  

Let $Lf$ and $Lg$ be the best affine approximations to $f$ and $g$ at
$(x_0, y_0)$.  Then
\begin{eqnarray*}
Lf(x,y) & = & 0\\
Lg(x,y) & = & 0
\end{eqnarray*}
is a system of {\it linear\/} equations.  If we are reasonably lucky, its 
coefficient matrix is invertible, so the system has a unique solution.  
Redefine $(x_0, y_0)$ to be this solution.  Iterate!  Here is a crude 
Maple implementation:

\begin{verbatim}
with(linalg):
f := x^3 + 3 * x^2 * y + y^3 + 2:
g := sin(x) + sin(y) + x * y + 5:
fx := diff(f,x):
fy := diff(f,y):
gx := diff(g,x):
gy := diff(g,y):
Lf := subs( x = x0, y = y0, f + fx * (X - x0) + fy * (Y - y0) ): 
Lg := subs( x = x0, y = y0, g + gx * (X - x0) + gy * (Y - y0) ):
x0 := 1.0:
y0 := 1.0:
for i from 1 to 10 do
     d := evalf( det( subs( x = x0, y = y0, array( [[fx,fy],[gx,gy]] ) ) ) ):
     if d = 0 then 
          print( `Try another initial guess.  `,
                 `The determinant is zero.` );
     fi;
     if abs(d) > 10^8 then
          print( `Try another initial guess.  Successive approximations `,
                 `appear to be diverging.` );
          break;
     fi;
     ans := fsolve( {Lf = 0, Lg = 0}, {X, Y} ):
     x0 := subs( ans, X ):
     y0 := subs( ans, Y ):
     print( x0, y0 );
od:
\end{verbatim}

\vspace{0.1in}

Here is the output:
\begin{verbatim}
8.64256  -11.6305
2.41747  -9.76917
1.17358  -6.69838
1.62284  -4.2193
2.15666  -2.26378
3.04366  -0.737431
3.31915  -1.18765
3.37967  -1.13957
3.37873  -1.14124
3.37873  -1.14124
\end{verbatim}
Observe the convergence!  What this gives us is a single approximate
solution $(3.37873, -1.14124)$ to the system of nonlinear equations
$$\setof{f(x,y) = 0, g(x,y) = 0}.$$%
Of course what you see printed out
are $10$ successive approximations, but it is the last (and best) 
approximation which is of interest.  If you want more than one approximate
solution, use another initial guess.

One thing you have to look out for is that an initial guess may not lead
to an approximate solution at all, no matter how many times you iterate.
If successive iterations yield bigger and bigger and bigger numbers, this
may be happening.  In the above code, the size of the determinant (of what?)
is checked.  Thus the determinant serves as an indicator of divergent
computation, and if it is too big, calculation terminates.

\vspace{0.1in}

\par\noindent Here are some Maple tools, used in the above example:

\vspace{0.1in}

\begin{tabular}{ll}
{\tt diff(f,x)} & compute the partial derivative of {\tt f} \WRT\ {\tt x}\\
{\tt evalf(x)} & give a numerical approximation to {\tt x}.
\end{tabular}

\vspace{0.1in}

\par\noindent Here is another very useful construction:
\begin{verbatim}
       sub( x = 1, y = 2, x^2 + y^2 + z^2 )
\end{verbatim}
which substitutes $1$ for $x$ and $2$ for $y$, yielding $5+z^2$.

\vspace{0.1in}

\par\noindent{\bf Exercises.}
Do this project in your group, turning in only one writeup for the whole group.
Problems 1 through 4 may be handwritten.  When you are done, each group member 
should turn in a peer evaluation sheet, as described on the first day sheet.

The first four exercises do not require use of the computer.  Be sure that your 
report includes the Maple code that you use, as well as judiciously selected 
portions of the output.

\vspace{0.1in}

\problemno{1} {\bf [4 points]}
Let \mp[[ f || \R^n || \R^k ]] and 
\mp[[ g || \R^k || \R^p ]] be mappings.  Then their composition is a mapping
\mp[[ g \circ f || \R^n || \R^p ]].
\begin{alphalist}
\item Use the alternate definition to show that if both $f$ and $g$ are
linear, then so is $g \circ f$.
\item Use the definition to show that if both $f$ and $g$ are
affine, then so is $g \circ f$.
\end{alphalist}

\vspace{0.1in}

\problemno{2} {\bf [3 points]}
Let \mp[[ f || \R^2 || \R ]].  Verify that the best affine 
approximation to $f$ at $(x_0,y_0)$ satisfies the {\it definition\/} 
of affine mapping.  What are $M$ and $c$?

\vspace{0.1in}

\problemno{3} {\bf [3 points]}
Let $f$ and $g$ be mappings from $\R^2$ to $\R$.  Fix
$(x_0,y_0) \in \R^2$.  If the system of linear equations
\begin{eqnarray*}
Lf(x,y) & = & 0\\
Lg(x,y) & = & 0
\end{eqnarray*}
is put in the form $Av = b$, where $v = \BRMAT{x\cr y}$, what are $A$ and $b$?

\vspace{0.1in}

\problemno{4} {\bf [3 points]} 
In the sample code shown for Newton's method, what 
is the purpose for checking to see if {\tt d = 0}?  

\vspace{0.1in}

\problemno{5}  {\bf [10 points]}
Find four approximate real solutions to the system
\begin{eqnarray*}
x^3 \sin(x) + y^2 \sin(y) + z^2 \sin(z) + x y z & = & 1\\
e^{1/(x^2 + y^2 + z^2 + 1)} - x - y - z & = & -1\\
x^2 y + y^2 z + \sin(x + y + z) & = & -9.
\end{eqnarray*}

Do this problem by copying with appropriate changes the given Maple
implementation of Newton's method.  You are to {\bf improve} the code,
as indicated in the following paragraph.

The given code always stops after $10$ iterations.  It is quite possible that
after $10$ iterations, $(x_0,y_0)$ is not at all close to a solution.  Or 
it could happen that $(x_0,y_0)$ is very close to a solution after less than
$10$ iterations.  {\bf Find a better stopping criterion, and implement it.}

{\bf Note.}  To get a number of real solutions, set the Maple
preferences so that old outputs are not always deleted.

\vspace{0.1in}

\problemno{6}  {\bf [12 points]}
Consider the system of equations:
\begin{eqnarray*}
-5x + 18x^2 + x^5 - 39y + 16xy + 20y^2 - 15xy^2 + 5xy^4 - 20\cos(xy) +
8e^{x+y} & = & 1\\
4x^5 + 13x^4 + 24x^2y + y^4 + x^3 y^3 + 4\cos(x+y) +
10e^{ 1 / (x^2 + y^2 + 0.2) }& = & 32.
\end{eqnarray*}
Use {\tt implicitplot} (see the worksheet {\tt newt.ms}) to create an
overlapping graph of the two equations.  Carefully adjust the ranges for
$x$ and $y$, if necessary generating more than one graph.  Use this information 
to guess the rough locations of {\bf all} common real solutions to the system.  
Then use Newton's method to get 
good estimates for all the common real solutions.  Also, find an approximate 
{\it complex\/} solution which is not a real solution.