Finding inscribed polygons on a curve.

We now take a look at the performance of Newton's method in a particular geometric application: finding inscribed polygons of various edge length ratios on a parametrized curve. We assume that we have a function

f[t_] := {(1 + 0.1*Cos[11*t])Cos[t],(1 + 0.1*Cos[11*t])Sin[t]}

Dist[t1_,t2_,t3_] := {(f[t2] - f[t1]).(f[t2] - f[t1]) - (f[t1] - f[0]).(f[t1] - f[0]),
(f[t3]  - f[t2]).(f[t3] - f[t2]) - (f[t2] - f[t1]).(f[t2] - f[t1]),
(f[0]     - f[t3]).(f[0] - f[t3]) - (f[t3] - f[t2]).(f[t3] - f[t2])}

f[0.2]
Dist[0.2,0.4, 0.6]
Dist[x,y,z]

We are now searching for a point where Dist[x,y,z] = {0,0,0}. This will be a point where the distances between consecutive points in our list are equal. To do so, we'll use Newton's method. We first need to compute the Jacobian of Dist.

MatrixForm[J=Outer[D,Dist[x,y,z],{x,y,z}]]

J /. {x->0.2,y->0.4,z->0.6}

We now set up the actual Newton iteration.

NStep[{s_,t_,u_}] := LinearSolve[J /. {x->s,y->t,z->u}, -Dist[s,t,u]] + {s,t,u};

NStep[{0.2,0.4,0.6}]
Column[NewtonSteps = NestList[NStep, {1.0,3.0,5.0}, 12]]

We now display the Newton steps graphically. To do so, we'll need a helper function MakePoints[{s,t,u}] which returns a list of the points on the curve which correspond to f[s], f[t], and f[u].

Bgplot := ParametricPlot[f[s],{s,
0,2*Pi},AxesOrigin->{0,0}];

MakePoints[{s_,t_,u_}] := {f[0],f[s],f[t],f[u],f[0]};

Manipulate[Show[Bgplot, ListLinePlot[MakePoints[NewtonSteps[[n]]], PlotStyle->{PointSize[Large]}, PlotMarkers->Automatic]],{n,1,12,1}]

Your book doesn't discuss this case, but Newton's method works just as well for underdetermined systems where there are fewer equations than variables. In this case, we need to find a least-squares solution to the matrix problem since in general the Jacobian matrix can't be nonsingular. 

Here's a worked example. 

Suppose first that we extend the previous code to find collections of 5 points in the plane which are at distance 1 from each other. We need to first create a new Dist function:

NewDist[X_,Y_,Z_,W_,U_] := {(Y-X).(Y-X) - 1, (Z-Y).(Z-Y) - 1, (W-Z).(W-Z) - 1, (U-W).(U-W) - 1, (X-U).(X-U) - 1}

Notice that X_, Y_, Z_, W_, and U_ are points in the plane, so there are 5 equations and 10 variables here.

MatrixForm[NewJ=Outer[D,NewDist[{x1,x2},{y1,y2},{z1,z2},{w1,w2},{u1,u2}],{x1,x2,y1,y2,z1,z2,w1,w2,u1,u2}]]

Again, the Newton step is found by solving the system Jh = -NewDist, but now we need to use Mathematica's LeastSquares solution. We try out Newton's method, starting with the collection of points (0,0), (1,0), (2,0), (2,1), (1,1). We will see that the method converges nicely.

NewNStep[{s1_,s2_,s3_,s4_,s5_,s6_,s7_,s8_,s9_,s10_}] := LeastSquares[NewJ /. {x1->s1,x2->s2,y1->s3, y2->s4, z1->s5, z2->s6, w1->s7,w2->s8,u1->s9,u2->s10}, -NewDist[{s1,s2},{s3,s4},{s5,s6},{s7,s8},{s9,s10}]] + {s1,s2,s3,s4,s5,s6,s7,s8,s9,s10};

NewNStep[{0.,1.,1.,1.,2.,1.,2.,2.,2.,2.}]

Column[NewNewtonSteps = NestList[NewNStep, {0.,1.,1.,1.,2.,1.,2.,2.,1.,2.}, 12]]

Again, displaying our results graphically requires a little bit of trickery in order to get the list of points to display right.

MakePoints[L_] := Join[Partition[L,2],{{L[[1]],L[[2]]}}];
MakePoints[NewNewtonSteps[[1]]]

Manipulate[ ListLinePlot[MakePoints[NewNewtonSteps[[n]]], PlotStyle->{PointSize[Large]}, PlotMarkers->Automatic, AxesOrigin->{0,0}],{n,1,12,1}]

Mathematica's Nest function allows us to just get the result of applying Newton's method for a certain number of steps (instead of listing all of the intermediate results).

Nest[NewNStep,{0.,1.,1.,1.,2.,1.,2.,2.,1.,2.}, 12]

Now comes the fun part! Suppose we want to drag this creature around by applying a force on the first vertex and the third vertex (say, up and down). We apply a little step in the direction of Velocity, then then apply Newton's method to correct the edgelengths after every step. This yields a fun visual effect.

h := 0.05;
Velocity[{x1_,x2_,x3_,x4_,x5_,x6_,x7_,x8_,x9_,x10_}] := {0,0,0,0,0.0,-0.01,0.0,0,0.0,0.01}; 
Velocity[{0.,1.,1.,1.,2.,1.,2.,2.,1.,2.}] 

(* Alternate Velocity function *)

(*Velocity[{x1_,x2_,x3_,x4_,x5_,x6_,x7_,x8_,x9_,x10_}] := 

Module[{c = 0.2*{x1+x3+x5+x7+x9,x2+x4+x6+x8+x10}},0.1*{x1 - c[[1]],x2 - c[[2]],x3 - c[[1]],x4 - c[[2]],x5 - c[[1]], x6 - c[[2]], x7 - c[[1]], x8 - c[[2]], x9 - c[[1]], x10 - c[[2]]}]; *)

(* *)


VArrows[{x1_,x2_,x3_,x4_,x5_,x6_,x7_,x8_,x9_,x10_}] := Module[ {V = Velocity[{x1,x2,x3,x4,x5,x6,x7,x8,x9,x10}]},

s1 = Sqrt[V[[1]]^2 + V[[2]]^2] + 0.01;
s2 = Sqrt[V[[3]]^2 + V[[4]]^2] + 0.01;
s3 = Sqrt[V[[5]]^2 + V[[6]]^2] + 0.01 ;
s4 = Sqrt[V[[7]]^2 + V[[8]]^2] + 0.01 ;
s5 = Sqrt[V[[9]]^2 + V[[10]]^2] + 0.01;

Graphics[{Arrow[{{x1,x2},{x1 + (1/s1)*V[[1]],x2 + (1/s1)*V[[2]]}}],
                  Arrow[{{x3,x4},{x3 + (1/s2)*V[[3]],x4 + (1/s2)*V[[4]]}}],
                  Arrow[{{x5,x6},{x5 + (1/s3)*V[[5]],x6 + (1/s3)*V[[6]]}}],
                  Arrow[{{x7,x8},{x7 + (1/s4)*V[[7]],x8 + (1/s4)*V[[8]]}}],
                  Arrow[{{x9,x10},{x9 + (1/s5)*V[[9]],x10 + (1/s5)*V[[10]]}}]
}]];

AnimationStep[{x1_,x2_,x3_,x4_,x5_,x6_,x7_,x8_,x9_,x10_}] := Nest[NewNStep,{x1,x2,x3,x4,x5,x6,x7,x8,x9,x10} + Velocity[{x1,x2,x3,x4,x5,x6,x7,x8,x9,x10}],20];

Column[AnimationSteps = NestList[AnimationStep,{0.,1.,1.,1.,2.,1.,1.5,2.,0.8,2.},360]];

Manipulate[ Show[VArrows[AnimationSteps[[n]]],Graphics[{Opacity[0.7],LightBlue,Polygon[MakePoints[AnimationSteps[[n]] ]]}],ListLinePlot[MakePoints[AnimationSteps[[n]]], PlotStyle->{PointSize[Large]}, PlotMarkers->Automatic, AxesOrigin->{0,0},AspectRatio->Automatic]],{n,2,360,1}]

A more mathematically interesting demonstration might try to maximize a function of the vertex positions subject to the constraint that the edgelengths are equal. To do this, we'd use the same animation step, but set the velocity to the gradient of the function that we're trying to maximize. An example is found in the commented out alternate velocity function above. 

However, that's a topic for later in the class!