# Curvature Application

I found these amazing animations at https://courses.engr.illinois.edu/tam212/avt.xhtml#avt which I think will be wonderfully useful when teaching multivariable calculus next semester. I've started to give Mathematica a try to emulate the first image.

d = 4;
r = 3;
l1 = d;
l2 = Pi*r;
l3 = d;
l4 = Pi*r;
totalLength = l1 + l2 + l3 + l4;
positionFunction[t_] := Module[{s, theta},
s = Mod[t, totalLength];
Which[
s < l1,
{-d/2 + s, 3},
s < l1 + l2,
theta = Pi/2 - (s - l1)/r;
{d/2 + r*Cos[theta], r Sin[theta]},
s < l1 + l2 + l3,
{d/2 - (s - (l1 + l2)), -3},
s < l1 + l2 + l3 + l4,
theta = 3 Pi/2 - (s - l1 - l2 - l3)/r;
{-d/2 + r*Cos[theta], r*Sin[theta]}
]
]
velocityFunction[t_] := Module[{x, theta},
s = Mod[t, totalLength];
Which[
s < l1,
{r, 0},
s < l1 + l2,
theta = Pi/2 - (s - l1)/r;
r {Sin[theta], -Cos[theta]},
s < l1 + l2 + l3,
{-r, 0},
s < l1 + l2 + l3 + l4,
theta = 3 Pi/2 - (s - l1 - l2 - l3)/r;
r {Sin[theta], -Cos[theta]}
]
]
accelerationFunction[t_] := Module[{x, theta},
s = Mod[t, totalLength];
Which[
s < l1,
{0, 0},
s < l1 + l2,
theta = Pi/2 - (s - l1)/r;
{-Cos[theta], -Sin[theta]},
s < l1 + l2 + l3,
{0, 0},
s < l1 + l2 + l3 + l4,
theta = 3 Pi/2 - (s - l1 - l2 - l3)/r;
{-Cos[theta], -Sin[theta]}
]
]
Manipulate[
Graphics[{
Line[{{-d/2, r}, {d/2, r}}],
Circle[{d/2, 0}, r, {Pi/2, -Pi/2}],
Line[{{d/2, -r}, {-d/2, -r}}],
Circle[{-d/2, 0}, r, {3 Pi/2, Pi/2}],
{Blue, Thickness[.01],
Arrow[{positionFunction[t],
positionFunction[t] + velocityFunction[t]}]},
{Magenta, Thickness[.01],
Arrow[{positionFunction[t],
positionFunction[t] + accelerationFunction[t]}]},
Red, PointSize[Large], Point[positionFunction[t]]
}, PlotRange -> {{-6, 6}, {-4, 4}}], {{t, 0.1}, 0, 100}]


Which produces this image.

I have a couple of questions.

1. I know you can slow the animation by clicking the minus icon at the end of the slide bar, then clicking the double down arrow several times until the speed is satisfactory. However, is there something you can put in the code to set the speed?

2. I still haven't been able to draw the bottom part of the first animation on https://courses.engr.illinois.edu/tam212/avt.xhtml#avt and note how it continues to move to the right (without everything jumping around). Don't think I've seen an example of this yet in Mathematica.

3. Finally, if anyone has better advice on how to reproduce this animation, it would surely be appreciated, because I think my students would benefit.

Thanks.

I used this image to understand Nikie's explanation of xTotal=15-x[15]/.sol:

However, is there something you can put in the code to set the speed?

Animate takes an option AnimationRate

I still haven't been able to draw the bottom part of the first animation on https://courses.engr.illinois.edu/tam212/avt.xhtml#avt and note how it continues to move to the right (without everything jumping around).

I think you can get as close as you want by playing with the options to Plot, ParametricPlot and Graphics. For example to prevent the coordinate system from "jumping around", give it an absolute PlotRange or combine the dynamic graphics with a static one using Show (which will then use the PlotRange and other options from the first graphic).

Finally, if anyone has better advice on how to reproduce this animation, it would surely be appreciated, because I think my students would benefit.

Not sure if it's better, but here's how I would do it:

tClothoid = 2/FresnelS[1];
(* define a piecewise linear function for the curvature *)
curvature[t_] :=
Piecewise[{
{0, t < 1},
{(t - 1)*4*(Pi/tClothoid^2), t < 1 + tClothoid/2},
{(4*Pi*(1 + tClothoid))/tClothoid^2 - t*4*(Pi/tClothoid^2), t < 1 + tClothoid},
{0, t < 10},
{1, t < 10 + Pi},
{0, True}}]

(* solve the differential equation *)
sol = First@
NDSolve[{
x[0] == 0, y[0] == 0, γ[0] == 0,
x'[t] == Cos[γ[t]],
y'[t] == Sin[γ[t]],
γ'[t] == curvature[t]},
{x, y, γ}, {t, 0, 20}];

tTotal = 15 - x[15] /. sol;
Animate[
Module[{tEnd = clock, currentPos, currentTangent, currentCurvature},
currentPos = {x[tEnd], y[tEnd]} /. sol;
currentTangent = {x'[tEnd], y'[tEnd]} /. sol;
currentCurvature = {x''[tEnd], y''[tEnd]} /. sol;
Column[{Show[
ParametricPlot[Evaluate[{x[t], y[t]} /. sol], {t, 0, tTotal},
ImageSize -> 500, PlotStyle -> Dotted, PlotRangePadding -> .5,
Epilog -> {Red, PointSize[Medium], Point[currentPos],
Arrow[{currentPos, currentPos + currentTangent}], Magenta,
If[Norm[currentCurvature] > 10^-5,
Arrow[{currentPos, currentPos + currentCurvature}]]}],
ParametricPlot[Evaluate[{x[t], y[t]} /. sol], {t, 0, tEnd},
ImageSize -> 500], Axes -> False],
Plot[curvature[t], {t, 0, tEnd},
PlotRange -> {{0, tTotal}, {0, 2}}, ImageSize -> 400,
Epilog -> {Red, PointSize[Medium],
Point[{tEnd, curvature[tEnd]}]}, Exclusions -> None,
AxesLabel -> {"t", "a"}, Ticks -> None]}]], {clock, 10^-5,
tTotal}, AnimationRate -> 0.1]


Can you point out an article that does a number of examples of the Fundamental Theorem of Plane Curves, Existence?

That's easy: Modern Differential Geometry of Curves and Surfaces with Mathematica. For example, here's the section about Clothoids and here's the section about assigning any curvature function to a curve (that's what I did above with NDSolve).

The main part I don't understand is how you came up with the curvature function in your opening steps.

I thought the question was about graphics, so I didn't include the derivation for curvature. It was a bit tedious, but not hard to understand:

First I want the clothoid to span 180° and I want a closed figure, so y[t] at the end of the clothoid should be 2, because that's the diameter of the circle curve on the left.

To make things as simple as possible for Solve, I just looked at the first 90° of the clothoid, where the curvature is simply linear:

curvature90[t_] := a*t;
sol90 = First@
DSolve[{x[0] == 0, y[0] == 0, γ[0] == 0,
x'[t] == Cos[γ[t]],
y'[t] == Sin[γ[t]], γ'[t] == curvature90[t]},
{x, y, γ}, t]


First, solve for the time when the integrated angle γ reaches 90°:

solT = Solve[{γ[t] == π/2} /. sol90, t]


$\left\{\left\{t\to -\frac{\sqrt{\pi }}{\sqrt{a}}\right\},\left\{t\to \frac{\sqrt{\pi }}{\sqrt{a}}\right\}\right\}$

Obviously, I want the second solution, where t is positive.

Then I can solve for a:

solA = Solve[{y[t] == 1} /. Last[solT] /. sol90, a]


and arrive at the time the clothoid should take for the first 90°:

Last[solT] /. solA


t -> 1/FresnelS[1]

so the time for the total 180° clothoid is tClothoid = 2/FresnelS[1];

I wanted to be able to play with different values for tClothoid, so I solved for a separately:

Solve[{γ[t] == π/2} /. sol90 /. t -> tClothoid/2, a]


a -> (4 π)/tClothoid^2

so that explains the (t - 1)*4*(Pi/tClothoid^2) in the piecewise expression.

For the decreasing curvature part, I simply wanted the two piecewise expressions from t<1+tClothoid/2 and t<1+tClothoid to meet at 1+tClothoid/2, so I solved for:

Solve[(t - 1)*4*(Pi/tClothoid^2) == b - t*4*(Pi/tClothoid^2) /.
t -> 1 + tClothoid/2, b]


and got b -> (4 π (1 + tClothoid))/tClothoid^2

So that's how I derived the piecewise function.

In retrospect, it probably would have been much easier to define it as:

tClothoid90 = 1/FresnelS[1];
curvature[t_] :=
UnitTriangle[(t - 4)/tClothoid90]/tClothoid90*π +
UnitBox[(t - 10)/π]


Finally, I want the figure closed, i.e. I want to repeat the animation at some point of time tTotal where x[tTotal] == x[0] == 0.

Since this is after the circular arc, x[t] is moving at unit speed in the x-direction, so x[t] == x[t0] + t - t0, where t0 is any value after the circle t0 > 10+Pi. Solve this for x[t] == 0 and you get t -> t0 - x[t0]. I chose t0 == 15 for no good reason (I think my original animation ended at t=15) and got

tTotal = 15 - x[15] /. sol;


(but I could have used 17.9 - x[17.9] /. sol or 10+Pi - x[10+Pi] /. sol and have gotten the same value of course)

• This is absolutely magnificent and something I would like to truly understand. The first difficulty is that most books or articles on differential geometry spend a lot of time on proofs of theory but do very few, if any, examples. I've done some searching online, but I haven't been able to find what I need to understand this answer. Can you point out an article that does a number of examples of the Fundamental Theorem of Plane Curves, Existence? – David Jul 2 '15 at 2:50
• The main part I don't understand is how you came up with the curvature function in your opening steps. Could you first explain tClothoid = 2/FresnelS[1]; What does FresnelS[1] return, and why divide 2 by the result? – David Jul 2 '15 at 2:54
• My next question would be a needed explanation of {(t - 1)*4*(Pi/tClothoid^2), t < 1 + tClothoid/2}, . What does this mean? It obviously gives the left side of the triangle in the lower figure and maybe the first half of the clothed in the upper figure, but why? Very hard for a rookie to understand. – David Jul 2 '15 at 2:57
• @nikie by chance have you written any books? if you wrote an "image/data processing/analysis in mathematica" book, i'd buy it (as long as it was not more than \$85). – amr Jul 2 '15 at 20:54
• @amr: Not yet. I'll think about it ;-) – Niki Estner Jul 3 '15 at 6:30

For Question 1: See the answer here. The number following AnimationRate should be the rate (per second) you want your parameter to increase.

For Question 2: You can do this via judicious use of GraphicsGrid, Max, and the PlotRange and AxesOrigin settings of Plot. If we define trange as the amount of time you want visible on your axes, then the following code does the trick:

trange = 30;
Manipulate[
GraphicsGrid[{{Graphics[{Line[{{-d/2, r}, {d/2, r}}],
Circle[{d/2, 0}, r, {Pi/2, -Pi/2}],
Line[{{d/2, -r}, {-d/2, -r}}],
Circle[{-d/2, 0}, r, {3 Pi/2, Pi/2}], {Blue, Thickness[.01],
Arrow[{positionFunction[t],
positionFunction[t] + velocityFunction[t]}]}, {Magenta,
Thickness[.01],
Arrow[{positionFunction[t],
positionFunction[t] + accelerationFunction[t]}]}, Red,
PointSize[Large], Point[positionFunction[t]]},
PlotRange -> {{-6, 6}, {-4, 4}}]},
{Plot[Norm[accelerationFunction[tau]], {tau, Max[0, t - trange], t},
PlotRange -> {{Max[0, t - trange], Max[t, trange]}, {0, 1.1}},
AxesOrigin -> {Max[0, t - trange], 0}, AspectRatio -> 1/3,

• Exclusions -> None ought to take care of the "breaks". – J. M. will be back soon Jul 1 '15 at 21:07