Newton’s method

A system of m equations in n unknowns can always be thought of as a single equation of the form
	F(x) = b
where F:R^n -> R^m is a function (with components F_i(x) equal to the left side of the ith equation), x is the unknown vector in R^n that we wish to solve for, and b is a vector in R^m with components b_i equal to the number of the right side of the ith equation.

For instance, if b1, b2 are given real numbers and we want to find the intersections between the vertical parabola y=x^2- b1 and the horizontal parabola x=y^2 - b2, then we must find the solutions (x,y) of the 2 x 2 system
	y - x^2 = b1
	x - y^2 = b2
which may be written F(x,y) = b where b=(b1,b2) and

F[{x_,y_}] = {y-x^2,x-y^2}

Try to solve this system exactly, and you'll likely end up trying to identify the roots of a degree four polynomial.  There's a known way to do this, but it's not pretty.   Take b = (-2,-3) for example.  The following picture clearly shows four solutions to the system:

ContourPlot[{y-x^2==-2,x-y^2==-3},{x,-5,5},{y,-5,5}]

And Mathematica will find them if you ask it to.

TableForm[{x,y}/.Simplify[Solve[F[{x,y}]=={-2,-3},{x,y}]]]

And this is about as easy as non-linear systems of equations get.  Typically, it's impossible to find exact solutions to non-linear systems. 

Nevertheless, when the function F is differentiable and there is a reasonably good initial guess x=x0 at a solution of F(x) = b, one can often use linear approximation to improve the guess repeatedly, thereby obtaining arbitrarily good approximations of an exact solution.  This is known as Newton's method.   Specificly, one replaces the function F in F(x) = b with its linear approximation at the initial guess x0 to get a new equation

F(x0) + Df(x0) (x-x0) = b

The great thing is that this is a linear system (as you can see by moving the F(x0) term to the other side of the equation), and we have a fairly straightforward way to solve linear systems.  To illustrate, let’s take the equation F(x,y) = (-2,-3) above and observe that if x0 = (2,2), then F(x0) = (-3,-3) is fairly close to our target vector.  We then employ Newton’s method as follows:

MatrixForm[DF[{x_,y_}] = D[F[{x,y}],{{x,y}}]]

x0 = {2.,2}
b={-2,-3}

Linapprox[{x_,y_}]= Simplify[F[x0]+DF[x0].({x,y}-x0)]

Our new guess is

x1={x,y}/.First[Solve[Linapprox[{x,y}]==b,{x,y}]]

Note that F(x1) is much closer to (-2,-3):

N[F[x1],25]

So already we have a lot of improvement.  However, the really great thing about Newton’s method is that you can repeat this step to improve your guess further:

Linapprox[{x_,y_}]= Simplify[F[x1]+DF[x1].({x,y}-x1)]
x1={x,y}/.First[Solve[Linapprox[{x,y}]==b,{x,y}]]
N[F[x1],25]

...and so on as often as you like

Nonlinear systems with more equations than unknowns

Both in theory and in practice, one often needs to solve (or at least understand solutions of) non-linear systems with more equations than unknowns.  Consider for instance the single equation f[x,y] = 0 where

f[{x_,y_}]=x^3+y^3-2x y

There are many points (x,y) that solve f(x,y)=0.  In fact the set of all such points, which includes (1,1) is a classical curve known as a `folium of Descartes’, which looks like this

Show[ContourPlot[f[{x,y}]==0,{x,-2,2},{y,-2,2}],Graphics[{PointSize->Large,Red,Point[{1,1}]}]]

The red point is (1,1).  From the picture it’s evident that if we vary either the x or the y coordinate of (1,1) a little bit then we can change the other coordinate a little bit to compensate and obtain another point on the curve.  For instance, we can hope that there is a point x near 1 such that (x,1..05) is still on the curve.  We can again use Newton’s method to help find this point.  But since we already have a value for y in mind, we might as well think of f as a function h(x) = f(x,1.05) as a function of x only

y0=.8
h[x_] =f[{x,y0}]

`

We want to solve h[x]=0 and again resort to linear approximation:

step = 0
x0 =1
Show[ContourPlot[f[{x,y}]==0,{x,-2,2},{y,-2,2}],Graphics[{PointSize->Large,Red,Point[{1,1}]}],
Graphics[{PointSize->Large,Blue,Point[{x0,y0}]}]]

step+=1
lapprox[x_] = h[x0]+h'[x0](x-x0)
x0=x/.First[Solve[lapprox[x]==0,x]]
h[x0]
Show[ContourPlot[f[{x,y}]==0,{x,-2,2},{y,-2,2}],Graphics[{PointSize->Large,Red,Point[{1,1}]}],
Graphics[{PointSize->Large,Blue,Point[{x0,y0}]}]]

The implicit function theorem says (in the present circumstance) that we can solve the equation f[x,y]=0 for x in terms of y near any point at with the partial derivative of f with respect to x is non-zero.  So let’s find the bad points:

fx[{x_,Y_}]=D[f[{x,y}],x]
badpts={x,y}/.Solve[{f[{x,y}]==0,fx[{x,y}]==0},{x,y}]
Show[ContourPlot[f[{x,y}]==0,{x,-2,2},{y,-2,2}],ListPlot[badpts,PlotStyle->{Red,PointSize[Large]}]]

The general idea is that if f[x1,...,xn] is a decent (i.e. C^1) scalar function of n variables and f[a] = b, then one can solve the equation f[x] = b for xj in terms of the other components of x near a provided the jth partial derivative of f at a is non-zero.  The 2 dimensional sphere in R^3 provides a nice illustration

f[{x_,y_,z_}]=x^2+y^2+z^2
grad[{x_,y_,z_}]=D[f[{x,y,z}],{{x,y,z}}]
plot1=ContourPlot3D[f[{x,y,z}]==1,{x,-2,2},{y,-2,2},{z,-2,2},AxesLabel->{x,y,z}];
plot2=ContourPlot3D[grad[{x,y,z}][[2]]==0,{x,-2,2},{y,-2,2},{z,-2,2},AxesLabel->{x,y,z}];
Show[plot1]

The torus in R^3 gives a more elaborate example.

R=2
r=1
f[{x_,y_,z_}]=(x^2+y^2+z^2+R^2-r^2)^2-4R^2 (x^2+y^2) 
grad[{x_,y_,z_}]=Simplify[D[f[{x,y,z}],{{x,y,z}}]]

plot1=ContourPlot3D[f[{x,y,z}]==0,{x,-3,3},{y,-3,3},{z,-3,3},AxesLabel->{x,y,z}];
plot2 =ContourPlot3D[grad[{x,y,z}][[1]]==0, {x,-3,3},{y,-3,3},{z,-3,3},AxesLabel->{x,y,z}];
Show[plot1]

Here is an example with two equations in three unknowns, an intersection of two cylinders.

f[{x_,y_,z_}]={x^2+y^2,(x-1)^2+z^2}
a={1,Sqrt[3],1}
b=f[a]
plot1=ContourPlot3D[{f[{x,y,z}][[1]]==b[[1]],f[{x,y,z}][[2]]==1},{x,-3,3},{y,-3,3},{z,-3,3}];
Show[plot1,Graphics3D[{PointSize->Large,Red,Point[a]}]]

Df[{x_,y_,z_}] =D[f[{x,y,z}],{{x,y,z}}]
vec=First[NullSpace[Df[a]]]
plot2=ParametricPlot3D[a+t vec,{t,-1,1},PlotRange->{{-3,3},{-3,3},{-3,3}},PlotStyle->{Thick}];
Show[plot1,plot2,Graphics3D[{PointSize[.03],Red,Point[a]}],AxesLabel->{x,y,z}]

MatrixForm[Df[{2,0,0}]]