Calculus IV with Maple
Copyright 2002, Dr. Jack Wagner
j.wagner@intelligentsearch.com

Lesson 5: Basic Differential Geometry of Curves

Topics: Construction of the tangent, normal and binormal (Frenet trihedron).  Curvature and torsion. The osculating circle. Construction of a moving Frenet trihedron  and a moving osculating circle. Projection of a spacecurve onto coordinate planes.

Maple commands introduced: Norm, Normalize, transform

Aside from the intrinsic interest in knowing the length of a particular curve, arc length is important as a parameter. The formulas of differential geometry take a particularly simple form when the curves involved are parameterized by arc length.  Because  ,  the computation of the derivative with respect to arc length, is the equivalent of normalizing .

In reading the code that  follows, observe the use of forward quotes (the upper left hand key on your keyboard).  The use of forward quotes creates a symbol.  Entries which otherwise would not be acceptable on the left side of an assignment operator may be used to create notation which resembles what we are accustomed to seeing.  

Example 5.1

Find the arc length for    between the points t=a and t=b.

>	restart:with(VectorCalculus): with(LinearAlgebra):

Warning, the assigned names <,> and <|> now have a global binding

Warning, these protected names have been redefined and unprotected: *, +, ., Vector, diff, int, limit, series

Warning, the names CrossProduct and DotProduct have been rebound

>	h := <t*cos(t), t*sin(t), t>;

>	dh[t] := diff(h, t);

Norm

>	ds[t] := Norm(dh[t], 2);

>	s := int(ds[t], t=a..b);

Here we see a difficulty with using the Norm  function.  Because of the absolute value functions (put in to accommodate complex variables), the int  function returns unevaluated.  We can get around this by telling Maple to assume t  is a real variable.

>	s := int(ds[t], t=a..b) assuming t::real;

Normalize

If you prefer to do this more directly without using Maple's Norm  function, you can enter  using the definition of the norm of dh[t], which is

>	ds[t] := sqrt(dh[t].dh[t]);

>	s := int(ds[t], t=a..b);

Maple's built-in function for arc length confirms the result.

ArcLength

>	ArcLength(h, t=a..b);

Example 5.2

We will now construct the Frenet trihedron to the same curve as in Example 5.1 at the point where t=5. That is, we will construct an orthogonal set of axes consisting of the unit tangent, the unit normal and the unit binormal to the curve at the point in question.  We need the unit tangent first.

>	restart: with(LinearAlgebra): with(VectorCalculus):with(plots):

Warning, the names CrossProduct and DotProduct have been rebound

Warning, the assigned names <,> and <|> now have a global binding

Warning, these protected names have been redefined and unprotected: *, +, ., Vector, diff, int, limit, series

Warning, the name changecoords has been redefined

Let Maple know we are dealing with real positive values of t.

>	assume( t > 0 ): interface(showassumed=0):

>	h := <t*cos(t), t*sin(t), t>;

>	dh[t] := diff(h, t);

>	ds[t] := Norm(dh[t], 2);

>	ds[t] := simplify( % );

The unit tangent.

>	dh[s] := dh[t]/ds[t];

>	T := eval(h+u*dh[s], t=5);

The Maple function Tangent  or TangentLine  would compute  

Tangent     

TangentLine

>	Tangent(h, t=5);

>	eval(h+u*dh[t], t=5);

>	TangentLine(h, t=5);

>	P1 := spacecurve(evalm(h), t=1..10, axes=boxed, color=blue):

>	P2 := spacecurve(evalm(T), u=-3..3, axes=boxed, color=red, labels=[x, y, z]):

>	display(P1, P2);

We have actually drawn the tangent three times unit length for ease of visualization.  

Now for the normal.  This is no different than the planar case except it is three dimensional. The derivative of  dh/ds which is a unit tangent vector, will be orthogonal to dh/ds. The unit normal is the normalized derivative.  .The  curvature of h at the point t, is defined by  =

>	dT[t] := simplify(diff(dh[s], t));

This non-normalized vector is produced by the Maple function.

PrincipalNormal

>	PrincipalNormal(h, t);

>	n := simplify(Normalize(dT[t], 2, inplace));

>	N := eval(h+u*n, t=5.0);

>	P3 := spacecurve(evalm(N), u=-3..3, axes=boxed, color=red, labels=[x, y, z]):

>	display(P1, P2, P3, scaling=constrained);

We easily calculate the curvature at t=5.

>	kappa(5) := evalf(Norm(subs(t=5, dT[t]/ds[t]), 2));

Alternatively,

Curvature

>	kappa(5) := eval(Curvature(h, t), t=5.0);

Now we need the binormal;  the unit vector orthogonal to both the tangent and normal and forming with them a right handed coordinate system i.e. the cross product. Note the order of the factors!

>	b := simplify(dh[s] &x n);

There is no need for normalization because   and N  and T are both unit orthogonal vectors.

>	B := evalf(subs(t=5, h+u*b));

Using the Maple built-in function Binormal:

Binormal

>	B := eval(h+u*Binormal(h, t), t=5.0);

>	P4 := spacecurve(evalm(B), u=-3..3, axes=boxed, color=red):

>	display(P1, P2, P3, P4);

Here we have the Frenet trihedron at t=5.  Referred to this set of coordinates the curve takes its simplest form in a neighborhood of t=5.  In fact, at a given point, p, the plane of the tangent and normal (the osculating plane) is the plane in which the curve lies.

Having exerted ourselves thusfar it is good to know that Maple has a built-in function to compute the Frenet frame.

TNBFrame

>

>	TNB := eval(TNBFrame(h, t), t=5.0);

>	p1 := eval(h, t=5.0);

>	Tngt := evalm(p1+u*TNB[1]);

>	Nml := evalm(p1+u*TNB[2]);

>	Bnml := evalm(p1+u*TNB[3]);

>	P5 := spacecurve({Tngt, Nml, Bnml}, u=-2..2):

>	display(P1, P5);

Example  5.3

We will find the equation of the curve in Example 6.2 in a neighborhood of the point .

In the adapted coordinates, = ''( a) + '''( a ) ... For small values of h, the second degree term will dominate the form of the curve.  This means that at any point, p= g(a), any curve referred to a local system of coordinates consisting of the Frenet trihedron will be almost a parabola in the osculating plane, in a neighborhood of p.  

We are going to omit the restart: that ordinarily precedes each new problem because we are going to need some of the results just obtained.

>	h := <t*cos(t), t*sin(t), t>;

We will work at the point where t=0.

>	p := eval(h, t=0);

>	dh[t] := diff(h, t);

>	dh[t=0] := subs(t=0, dh[t]);

The tangent lies in the x-z plane and makes an angle of    with the x axis.

>	dh[tt] := diff(h, t$2);

>	dh[tt=0] := subs(t=0, dh[tt]);

The normal is parallel to the y axis.  The osculating plane makes an angle of  with the x-y plane. First we calculate the first three terms of the Taylor expansion.  Remember, the first term is zero.

>	g := dh[t=0]*t+dh[tt=0]*t^2/2;

Now we will rotate the equation for the curve by  so that the osculating plane will be parallel to the x-y plane, and the Taylor expansion approximation, "g",   will also be rotated by the same amount.  This has the effect of making the Frenet trihedron the coordinate axis vectors.  We will need the rotation matrix about the y-axis.

>	R := <<cos(theta)|0|-sin(theta)>, <0|1|0>, <sin(theta)|0|cos(theta)>>;

>	R := subs(theta=-Pi/4, R);

The reader familiar with the elements of linear algebra will recognize this as the rotation matrix which, when applied to a vector,  will rotate it by    around the y axis. Those who do not have the pleasure of this familiarity will have to take it on faith.

>	ht := <<h[1]>, <h[2]>, <h[3]>>;

>	H := R.ht;

>	g := <<g[1]>, <g[2]>, <g[3]>>;

>	g2 := R.g;

>	S1 := spacecurve(evalm(H), t=-.05..0.05, color=blue):

>	S2 := spacecurve(evalm(g2), t=-.05..0.05, color=red):

>	display(S1, S2, labels=[x, y, z], axes=normal);

This illustrates the fact that in a local coordinate system, there is a neighborhood of the origin in which the curve is approximated by a parabola in the osculating plane.  

In addition to the curvature, there is another important invariant of the curve, the torsion, which, as its name implies, measures the twisting of the curve in space.  Torsion is computed from the derivative of the binormal vector.  It is, in fact, parallel to, but oppositely oriented from, the normal.  

That is   , where the torsion, , is a differentiable function. The torsion is computed as:

>	db[s] := simplify(diff(b, t)/ds[t]);

>	tau(5) := -eval(db[s].n, t=5.0);

Torsion

>	eval(Torsion(h, t), t=5.0);

More generally, we may refer a curve to axes formed by the Frenet frame at any given point. Consider the Taylor expansion, to third order,  of a function f .   .

In the basis formed by the Frenet frame, ,    ,  ,

>	restart: with(LinearAlgebra): with(VectorCalculus):with(plots):

Warning, the names CrossProduct and DotProduct have been rebound

Warning, the assigned names <,> and <|> now have a global binding

Warning, these protected names have been redefined and unprotected: *, +, ., Vector, diff, int, limit, series

Warning, the name changecoords has been redefined

>	f[s] := <1, 0, 0>;

>	f[ss] := <0, kappa, 0>;

>	f[sss] := <-kappa^2, d_kappa, kappa*tau>;

>	F := s*f[s]+s^2/2*f[ss]+s^3/6*f[sss];

>	interface(showassumed=0):

>	assume(t>0):

>	h := <t*cos(t)|t*sin(t)|t>;

>	dh[t] := diff(h, t);

>	kappa := Curvature(h, t):

>	tau := Torsion(h, t):

>	d_kappa := diff(kappa, t)/Norm(dh[t], 2):

>	kappa := eval(kappa, t=5.0);

>	tau := eval(tau, t=5.0);

>	d_kappa := eval(d_kappa, t=5.0);

>

F := F;

>	spacecurve(evalm(F), s=-5..5, axes=normal, labels=[t, n, b], color=blue, title="The curve h(t) plotted on the Frenet frame at t=5");

Having gone this far, there is one more interesting construction we would like to make.  At each point on a curve, one can draw a circle, lying in the plane of the tangent and the normal, the osculating plane, which has the same curvature as the curve, the so called circle of curvature. The radius of this circle, the radius of curvature,  is computed as  .

Example  5.4

Find and plot the circle of curvature at t=5.
Recollect that the polar representation of a circle centered at a point (0, p) with radius r, is given by    .  Referred to the Frenet  coordinate system in the osculating plane, this translates to     , for a circle of radius , centered r units along the normal.  Therefore, we plot:

We multiply by a factor of three to improve visibility

>	rho := 3*1/kappa(5);

>	C := <rho*cos(theta)| rho+rho*sin(theta)|0>;

>	h_5 := subs(t=5, h);

>	CC := C+h_5;

>	P1 := spacecurve(evalm(h), t=0..10):

>	P5 := spacecurve(evalm(CC), theta=0..2*Pi, axes=boxed, color=red):

>	display(P1, P5);

By rotating this plot you will be able to convince yourself that the circle matches the curvature of the curve at the point of tangency.

What we would really like to see is curvature and torsion at work.  What do they really look like?  After all, they are described by derivatives which imply motion.  We need to put all of this together, compute a sequence of plots of the Frenet trihedron at different points on the curve and then display these plots in sequence to get a "moving picture".

The following code accomplishes this.  To make animations work, you must right click on the plot, select animation from the menu and click on "play".  The animation may be toggled between single cycle and continuous and may be slowed down or speeded up from this submenu.  Be forewarned!  The sequences will take some 10 to 30 seconds each to be processed.

>	restart: with(LinearAlgebra): with(VectorCalculus): with(plots):

Warning, the names CrossProduct and DotProduct have been rebound

Warning, the assigned names <,> and <|> now have a global binding

Warning, these protected names have been redefined and unprotected: *, +, ., Vector, diff, int, limit, series

Warning, the name changecoords has been redefined

>	interface(showassumed=0):

>	assume(t>0):

>	h := <t*cos(t), t*sin(t), t>;

>	TNB := TNBFrame(h, t):

>	Tn := j->s*subs(t=j, TNB[1]):

>	Nn := j->s*subs(t=j, TNB[2]):

>	Bn := j->s*subs(t=j, TNB[3]):

>	H := j->subs(t=j, h);

>	TT := evalm(H(j/10)+Tn(j/10)):

>	S1 := seq(spacecurve(TT, s=-2..2, color=red), j=1..100):

>	NN := evalm(H(j/10)+Nn(j/10)):

>	S2 := seq(spacecurve(NN, s=-2..2, color=red), j=1..100):

>	BB := evalm(H(j/10)+Bn(j/10)):

>	S3 := seq(spacecurve(BB, s=-2..2, color=red), j=1..100):

>	d1 := display(S1, insequence=true):  d2 := display(S2, insequence=true):  d3 := display(S3, insequence=true):  d4 := spacecurve([h[1], h[2], h[3]], t=0..10, color=blue): display(d1, d2, d3, d4, scaling=constrained, orientation=[20,75]);

In some ways even more fascinating is the moving circle of curvature.  Its twists, turns and changing size clearly demonstrate the meaning of torsion and curvature.  Remember, it lives in the osculating plane which also contains the derivative of the binormal.  Here is the code to produce the moving circle of curvature.

>	restart: with(LinearAlgebra): with(VectorCalculus):with(plots):

Warning, the names CrossProduct and DotProduct have been rebound

Warning, the assigned names <,> and <|> now have a global binding

Warning, these protected names have been redefined and unprotected: *, +, ., Vector, diff, int, limit, series

Warning, the name changecoords has been redefined

>	assume(t>0):

>	h := <t*cos(t), 2*t*sin(t), t>:

>	dh[t] := diff(h, t):

>	ds[t] := Norm(dh[t], 2):

>	dh[s] := dh[t]/ds[t]:

>	dT := diff(dh[s], t):

>	n := simplify(dT/ds[t]):

>	rho := j->evalf(1/subs(t=j, Norm(n, 2))):

>	H := i->subs(t=i, h):

>	T := p->subs(t=p, dh[s]):

>	N := q->subs(t=q, n/Norm(n, 2)):

>	Z := H(k/10)+rho(k/10)*(N(k/10)+(cos(theta)*T(k/10)+sin(theta)*N(k/10))):

>	Circ := k->evalm(Z):

>	S := seq(spacecurve(Circ(k), theta=0..2*Pi), k=0..100):

>	P1 := display(S, insequence=true, color=blue):

>	P2 := spacecurve([h[1], h[2], h[3]], t=0..10, color=red):

>	display(P1, P2, axes=framed);

     In both these code sequences, we started off with assume(t>0). This greatly simplifies the form of the results obtained and cuts down on subsequent computation time.
It is sometimes of interest to project a spacecurve (or a surface) onto a coordinate plane or onto a plane parallel to a coordinate plane. This can be accomplished with the transform function in the plottools package.

transform

        To use the transform function which operates on plot structures it is first necessary to define the plot structure and then specify the parameters  to be transformed.

>	restart: with(plottools):  with(plots):  with(LinearAlgebra): with(VectorCalculus):

Warning, the names arrow and changecoords have been redefined

Warning, the names CrossProduct and DotProduct have been rebound

Warning, the assigned names <,> and <|> now have a global binding

Warning, these protected names have been redefined and unprotected: *, +, ., Vector, diff, int, limit, series

>	h := <t*cos(t)|2*t*sin(t)|t>:

>	c1 := spacecurve(evalm(h), t=1..10, axes=boxed, color=red):

>	T1 := transform((x, y, z)->[x, y, 0]):

          This is a tricky syntax.  The parameters to be transformed are in round brackets, and the parameters to which they are to be transformed are in square brackets.  Of course, we are not limited to coordinate projections.  A variety of parameter changes could be implemented, including, in particular, coordinate system transformations.  

>	c2 := spacecurve(evalm(h), t=1..10, axes=boxed, color=blue):

        We needed a second definition of the curve in a different color. One of the plot structures, c1, now becomes an argument to g, and is displayed with the second plot structure, c2.  The original curve and its projection onto the x-y plane are then displayed.

>	display(T1(c1), c2);

To project onto the x-z plane we enter:

>	T2 := transform((x, y, z)->[x, 0, z]):

>	display(T2(c1), c2);

>

Practice

2. For each of the curves in prob.#1, compute the curvature at the point indicated and plot the osculating circle and Frenet trihedron.

3. Select one of the curves in prob#1 and create an animation for;
       a) The moving Frenet trihedron
       b) The moving osculating circle.

4. Project the curves in Prob#1 onto  the coordinate planes and plot each of these projections together with the curve.