Newton’s Method in N Dimensions.

Step 1. Deriving the iteration.

Suppose we have a function (OverVector[F]): (ℝ^n)-> (ℝ^n), and we want to find points (OverVector[x]) where 

(OverVector[F])((OverVector[x])) = (OverVector[0])

Following the general idea of Newton’s method, we’ll observe that if we’re at a point (Subscript[(OverVector[x]), i]) and we take a step (OverVector[h]) we expect to end up at a point where 

(OverVector[F])((Subscript[(OverVector[x]), i])+ (OverVector[h])) ≃ (OverVector[F])((Subscript[(OverVector[x]), i])) + D (OverVector[F])((Subscript[(OverVector[x]), i]))(OverVector[h]) 

where D (OverVector[F]) is the Jacobian matrix of first partials of (OverVector[F]). That means that if we make 

(OverVector[h])= -(((D (Subscript[(OverVector[F]), (Subscript[(OverVector[x]), i])])))^(-1))((OverVector[F])((Subscript[(OverVector[x]), i])))

then

(OverVector[F])((Subscript[(OverVector[x]), i])+ (OverVector[h])) ≃ (OverVector[F])((Subscript[(OverVector[x]), i])) + D (OverVector[F])((Subscript[(OverVector[x]), i]))(OverVector[h])
   = (OverVector[F])((Subscript[(OverVector[x]), i])) + D (OverVector[F])((Subscript[(OverVector[x]), i]))(-(((D (Subscript[(OverVector[F]), (Subscript[(OverVector[x]), i])])))^(-1))((OverVector[F])((Subscript[(OverVector[x]), i]))))
   = (OverVector[F])((Subscript[(OverVector[x]), i])) - (OverVector[F])((Subscript[(OverVector[x]), i])) = (OverVector[0]). 

If we call (Subscript[(OverVector[x]), i ])+ (OverVector[h])= (Subscript[(OverVector[x]), i+1]), we can rewrite the Newton iteration as

(Subscript[(OverVector[x]), i+1 ])=(Subscript[(OverVector[x]), i])-(((D (Subscript[(OverVector[F]), (Subscript[(OverVector[x]), i])])))^(-1))((OverVector[F])((Subscript[(OverVector[x]), i]))) 

which really is the vector analogue of the scalar Newton iteration

(Subscript[x, i+1])= (Subscript[x, i ])- ((f((Subscript[x, i])))/(f'((Subscript[x, i]))))

especially if we (deliberately) confuse the notation for “-1st power” and “inverse” and write the scalar Newton iteration as

(Subscript[x, i+1 ])= (Subscript[x, i])- (((f'((Subscript[x, i]))))^(-1))f((Subscript[x, i])).

Part II. A worked example (with algebra).

Let’s write this out for a particular example. The complex polynomial 

f(z) = (z^3)- 1

can be thought of as a map (ℝ^2) -> (ℝ^2). We can write this out by taking:

cplxForm = Expand[(a+b I)^3 - 1]

Separating this into Real and Imaginary parts gives us the coordinate functions of (OverVector[F]). We can get the real and imaginary parts of the expression with ReIm, but be careful-- we have to tell Mathematica that the variables are real by wrapping the whole thing in those assumptions and simplifying

Assuming[{a ∈ Reals, b ∈ Reals},FullSimplify[ReIm[cplxForm]]]

We can finally define the function!

F[{a_,b_}] := {-1+(a^3)-3 a (b^2),3 (a^2) b-(b^3)};

Now we compute the inverse Jacobian. To do that, we’ll use the D function.

(DF = D[F[{a,b}],{{a,b}}]) // MatrixForm

We can invert the matrix (symbolically) using Inverse:

Inverse[DF] // MatrixForm

It’s easy to check that this worked:

Simplify[Inverse[DF].DF] // MatrixForm

Now we’re finally ready to construct the Newton iteration.

G[{a_,b_}] := {a,b} - (({{((3 (a^2)-3 (b^2))/(9 (a^4)+18 (a^2) (b^2)+9 (b^4))), ((6 a b)/(9 (a^4)+18 (a^2) (b^2)+9 (b^4)))}, {-((6 a b)/(9 (a^4)+18 (a^2) (b^2)+9 (b^4))), ((3 (a^2)-3 (b^2))/(9 (a^4)+18 (a^2) (b^2)+9 (b^4)))}})).F[{a,b}];

We can now check that this works for the solutions we know about

Solutions = {{1,0},{-1/2,(Sqrt[3])/2},{-1/2,-(Sqrt[3])/2}}

G /@ Solutions

Now we’ll try to use the iteration.

NestList[G,{2.0,3.0},17] // TableForm

Plotting the solution shows you that this iteration hops around quite unpredictably before settling down to converge to the solution:

SolutionsGraphic= Graphics[{{Red,PointSize[0.025],Point[Solutions[[1]]]},{Green,PointSize[0.025],Point[Solutions[[2]]]},
{Blue,PointSize[0.025],Point[Solutions[[3]]]}}];

Show[SolutionsGraphic,
Graphics[
Line[NestList[G,{2.0,3.0},22]]
],Axes->True,PlotRange->{{-6,6},{-6,6}}]

The dependence on the initial position is quite sensitive. To see it, we’ll write a little function which returns the color corresponding to the solution we get. We can iterate a nest “until it converges” with FixedPoint:

FixedPoint[G,{2.0,3.0}]

SolutionColor[{a_,b_}] := 
First[Nearest[Solutions,FixedPoint[G,{a,b}]]] /. {{1,0}->Red,{-1/2,(Sqrt[3])/2}-> Green, {-1/2,-(Sqrt[3])/2}-> Blue};

SolutionColor[{2.0,3.0}]

A quick and dirty way to put together boxes of these colors is to use ArrayPlot:

Block[{σ = 200},
ArrayPlot[
ParallelMap[SolutionColor,Table[{x+1.0/(2σ),y+1.0/(2σ)},{y,-1,1,1/σ},{x,-1,1,1/σ}],{2}]
]
]

If you have some cycles to burn, you can make a higher resolution plot by changing the resolution σ to a larger number.

If you start Newton’s method in the middle of a solid colored region, it will converge to the corresponding root. But what about the boundaries between the regions? These are actually a weird, lacy set of points where Newton’s method will not converge to any of the three roots. Those points are called a Julia set, and there’s a built-in way to make nice plots of them.

If we write down the Newton iteration explicitly as a rational function, we can feed it to JuliaSetPlot and get:

JuliaSetPlot[z - ((z^3)-1)/(3(z^2)),z,PlotRange->{{-1,1},{-1,1}},PlotTheme->"Web"]

A more applied example.

Here’s a more motivated example. Suppose that we have a robot arm in the plane with two rotary joints. We’d like to solve the “inverse kinematics” problem of finding joint angles which result in the end of the arm being at a particular target location ((Subscript[t, x]),(Subscript[t, y])) in the plane. 

We know that the end position of the arm is given by

ArmPos[{θ1_,θ2_}] := {Cos[θ1],Sin[θ1]} + {Cos[θ2],Sin[θ2]};

So we’re trying to solve the problem ArmError = (OverVector[0]) where

ArmError[{θ1_,θ2_},{tx_,ty_}] := ArmPos[{θ1,θ2}] - {tx,ty};

D[ArmError[{θ1,θ2},{tx,ty}],{{θ1,θ2}}]

Now in theory this is still simple enough to invert symbolically. But in practice, you may not be given the Jacobian of the function as a formula, just as piece of code that you can call. To handle this, we return to the basic equation

(OverVector[h])= -(((D (Subscript[(OverVector[F]), (Subscript[(OverVector[x]), i])])))^(-1))((OverVector[F])((Subscript[(OverVector[x]), i])))

and realize that if we apply D (Subscript[(OverVector[F]), (Subscript[(OverVector[x]), i])]) to both sides of the equation, we get the linear system of equations:

(D (Subscript[(OverVector[F]), (Subscript[(OverVector[x]), i])])) (OverVector[h]) = -(OverVector[F])((Subscript[(OverVector[x]), i])) 

This means that given

ArmErrorJacobian[{θ1_,θ2_},{tx_,ty_}] := {{-Sin[θ1],-Sin[θ2]},{Cos[θ1],Cos[θ2]}};

we can write the Newton iteration as

ArmNewtonIteration[{θ1_,θ2_},{tx_,ty_}] := {θ1,θ2} + LinearSolve[ArmErrorJacobian[{θ1,θ2},{tx,ty}],-ArmError[{θ1,θ2},{tx,ty}]];

Just to check it out, we try reaching for a spot that ought to be well within range: {1.2,1.2}.

Nest[ArmNewtonIteration[#,{1.2,1.2}]&,{0.0,1.0},19]

ArmPos[%]

Now the beauty of Newton iteration is that it’s very good at adjusting a solution as the problem changes continuously, after all, the current solution to the problem is an excellent starting point for Newton. So we can first define:

ArmUpdate[{θ1_,θ2_},{tx_,ty_}] := 
Nest[ArmNewtonIteration[#,{tx,ty}]&,{θ1,θ2},20];

This converges in a few thousandths of a second on my laptop

ArmUpdate[{0.0,1.0},{1.2,1.2}] // Timing

Now we’ll make a little demonstration using this principle. First, we’ll draw the robot arm from the given angles:

ArmGraphic[{θ1_,θ2_}] := Graphics[{FaceForm[Gray],EdgeForm[Black],
Rotate[Rectangle[{-0.05,-0.05},{1.05,0.05},RoundingRadius->0.05],θ1,{0,0}],
Translate[Rotate[Rectangle[{-0.05,-0.05},{1.05,0.05},RoundingRadius->0.05],θ2,{0,0}],{Cos[θ1],Sin[θ1]}],
{PointSize[Large],Black,Point[ArmPos[{θ1,θ2}]],Point[{Cos[θ1],Sin[θ1]}],Point[{0,0}]},{Dotted,Red,Circle[{0,0},2]}},Axes->True,PlotRange->{{-2,2},{-2,2}}];

Here’s our initial solution:

ArmUpdate[{0.0,1.0},{1.2,1.2}]

ArmGraphic[ArmUpdate[{0.0,1.0},{1.2,1.2}]]

Now we can update this as we change the location of the point. We’ll start by setting some base angles.

θold = {0.0,1.0}

Now we can write a Manipulate which updates the arm angles every time we move the point:

Manipulate[
θnew = ArmUpdate[θold,pt];
θold= θnew;
ArmGraphic[θnew],
{{pt,{1.2,1.2}},Locator}]