Newton’s method for non-linear systems

Here I implement Newton’s method to solve a non-linear system of two equations in two unknowns.  What’s needed to get started is
 
	a mapping 		 		F:R^2->R^2
	an `initial’ point in the domain:		guess0 
	a ‘target’ point 	in the range		targ

If F is differentiable, then one can define the `Newton map’

	N(p) = p + JF(p)^(-1)( targ - F(p) )

As one can check, a fixed point of N(p) is a solution of the system 

	F(p) = targ

And in fact, if F is, say, C^2 and the target point is close enough to F(guess0), then one can show that the sequence 

	guess0, N(guess0), N^2(guess0), ...

converges like gangbusters to a fixed point.

Below is the example I did in class.  Note that I enclose each command in `MatrixForm’ just to format the output more nicely.  You can easily change F, guess0, targ, to deal with other systems.  It requires a bit more modification to change the number of equations/unknowns to something other than 2.

MatrixForm[F[{x_,y_}] = {x^2-y,y^2-x}]
MatrixForm[JF[{x_,y_}]=Transpose[{D[F[{x,y}],x],D[F[{x,y}],y]}]]
MatrixForm[guess0={2,3}]
MatrixForm[F[guess0]]
MatrixForm[targ = {1/2,8}]

...and here we set up the Newton map.  For the sake of comparison, I also set up the ‘relaxed’ Newton map used in the proof of the inverse function theorem.

MatrixForm[Newt[{x_,y_}]=Simplify[{x,y}-Inverse[JF[{x,y}]].(F[{x,y}]-targ)]]
MatrixForm[RelaxedNewt[{x_,y_}]=Simplify[{x,y}-Inverse[JF[guess0]].(F[{x,y}]-targ)]]

Now I’m ready to start trying to solve my system.  The variable ‘n’ here keeps track of the number of iterates, and guess keeps track of our current approximate solution of the system.

Initially, we have

n=0
guess = guess0
F[guess]

Then we can execute the following block as often as necessary to improve our initial approximate solution.   (Enclosing a command in N[...,20] tells Mathematica to give the output in decimal form to 20 places.)

n=n+1
MatrixForm[N[guess =Newt[guess],20]]
MatrixForm[N[F[guess],20]]

For the sake of honesty, I should point out that Mathematica provides the ‘Nest’ and `NestList’ commands to make iterating functions easier.  

Here I use `NestList’ to compare Newton’s method with ‘relaxed’ Newton’s method.  It’s often the case in math that the best method to use in practice is rather different than the easiest method to deal with in theory.

N[TableForm[NestList[Newt,guess0,10]],20]

N[TableForm[NestList[RelaxedNewt,guess0,20]],20]