Diff_Geom_Curves.mws
Calculus with Maple
by Dr. Jack Wagner, Lehman College, City University of New York
Copyright 2001
Basic Differential Geometry of Curves
A
r
c length. 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.
map, norm, normalize, transform
We begin by looking at curves embedded in
. A curve is a map,
:
A segment of the real line is bent and twisted. It is those bends and twists that we will examine. Curves are ordinarily in the form:
. This is a frequently also written as
. Both these forms suggest that the curve is traced out by the tip of a vector whose coordinates are the three functions in brackets. The derivative is then simply found by differentiating each of the three functions. A curve which is continuous and has continuous non-vanishing derivatives is known as a regular curve. We will be working with regular curves. When a curve (or a surface) is not regular the usual rules of calculus fail. From time to time we will have occasion to point out examples. First, note that a regular curve is rectifiable. Its length can be measured by a limiting process on inscribed polygons. In three dimensions the length of a curve from t=a to t=b is:
=
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.
Ex. 6.1
Find the arc length for x(t)=[t*cos(t),t*sin(t),t] between the points t=a and t=b.
> |
restart: with(plots): with(linalg):
|
Warning, new definition for norm
Warning, new definition for trace
> |
h:=vector([t*cos(t),t*sin(t),t]);
|
> |
`dh/dt`:=map(diff,h,t);
|
map
norm
> |
`ds/dt`:=norm(`dh/dt`,2);
|
> |
s:=int(`ds/dt`,t=a..b);
|
Here we see the problem with using the norm function. Because of the absolute value functions (put in to accommodate complex variables), the int function returns unevaluated. The normalize function has the same problem.
normalize
Both these functions are useful in dealing with numerical calculations.
Back to basics!
> |
`ds/dt`:=sqrt(innerprod(`dh/dt`,`dh/dt`));
|
> |
s:=int(`ds/dt`,t=a..b);
|
Ex. 6.2
> |
We will now construct the Frenet trihedron to the same curve as in Ex. 6.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(plots): with(linalg):
|
Warning, new definition for norm
Warning, new definition for trace
> |
assume(t>0):#Let Maple know we are dealing with real positive values of t.
|
> |
h:=vector([t*cos(t),t*sin(t),t]);
`dh/dt`:=map(diff,h,t);
`ds/dt`:=sqrt(innerprod(`dh/dt`,`dh/dt`));
`dh/ds`:=evalm(eval(`dh/dt`)/eval(`ds/dt`)); #The unit tangent.
T:=evalm(subs(t=5,eval(h))+u*subs(t=5,eval(`dh/ds`)));
|
The eval is critical for reading the various symbols as vectors or lists.
> |
P1:=spacecurve(h,t=1..10,axes=boxed,color=blue):
|
> |
P2:=spacecurve(T,u=-3..3,axes=boxed,color=red,labels=[x,y,z]):
|
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/dt`:=simplify(map(diff,
eval(`dh/ds`),t));
|
> |
n:=simplify(evalm(eval(`dT/dt`)/sqrt(innerprod (eval(`dT/dt`),eval(`dT/dt`)))));
|
> |
N:=evalm(subs(t=5,eval(h))+u*subs(t=5,eval(n)));
|
> |
P3:=spacecurve(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,eval(`dT/dt`)/eval(`ds/dt`))));
|
In order to perform a substitution in a list or vector, we must first apply the eval operator.
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.
> |
b:=crossprod(eval(n),eval(`dh/ds`));
|
There is no need for normalization because
and N and T are both unit orthogonal vectors.
> |
B:=evalm(evalf(subs(t=5,eval(h))+u*subs(t=5,eval(b))));
|
> |
P4:=spacecurve(B,u=-3..3,axes=boxed,color=red):
|
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.
Ex. 6.3
We will find the equation of the curve in Ex. 6.2 in a neighborhood of the point
.
By a suitable translation of the origin to p, we can always arrange things so
In the coordinates represented by the Frenet trihedron, the curve is tangent to the "x" axis:
'(0)=0. Recollect the Taylor expansion for functions of a single variable:
=
+
'(
a)h+
''(
a
)
+...
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:=evalf(subs(t=0,eval(h)));
|
> |
`Dh[t]`:=map(diff,h,t);
|
> |
`Dh[t](0)`:=evalf(subs(t=0,eval(`Dh[t]`)));
|
The tangent lies in the x-z plane and makes an angle of
with the x axis.
> |
`Dh[t,t]`:=map(diff,h,t$2);
|
> |
`Dh[t,t](0)`:=evalf(subs(t=0,eval(`Dh[t,t]`)));
|
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:=evalm(`Dh[t](0)`*t+`Dh[t,t](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:=matrix(3,3,[cos(theta),0,-sin(theta),0,1,0,sin(theta),0,cos(theta)]);
|
> |
R:=evalf(subs(theta=-Pi/4,eval(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.
> |
S1:=spacecurve(H,t=-.05..0.05,color=blue):
|
> |
S2:=spacecurve(g2,t=-.05..0.05,color=red):
|
> |
display(S1,S2,labels=[x,y,z],axes=normal);
|
This illustrates the fact that in an adapted 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.
Maple computes this quantity easily as the 2-norm (the ordinary Euclidean norm) of the derivative of b.
> |
`db/dt`:=simplify(map(diff,eval(b),t));
|
> |
tau(5):=evalf(norm(subs(t=5,eval(`db/dt`))),2);
|
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
.
Ex. 6.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:
> |
rho:=3*1/kappa(5); # We multiply by a factor of three to improve //visibility
|
> |
C:=[rho*cos(theta), rho+rho*sin(theta),0];
|
> |
P5:=spacecurve(evalm(subs(t=5,eval(h))+C),theta=0..2*Pi,axes=boxed,color=red):
|
> |
display(P1,P2,P3,P4,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(plots): with(linalg):
|
Warning, new definition for norm
Warning, new definition for trace
> |
assume(t>0):
h:=[t*cos(t),t*sin(t),t]:
`dh/dt`:=map(diff,h,t):
T:=evalm(`dh/dt`/sqrt(innerprod(`dh/dt`,`dh/dt`))):
`dT/dt`:=map(diff,T,t):
N:=eval(`dT/dt`)/sqrt(innerprod(`dT/dt`,`dT/dt`)):
B:=crossprod(T,N):
Tn:=j->evalm(s*(subs(t=j,eval(T)))):
Nn:=j->evalm(s*(subs(t=j,eval(N)))):
Bn:=j->evalm(s*(subs(t=j,eval(B)))):
H:=j->(subs(t=j,eval(h))):
S1:=seq(spacecurve(evalm(H(j/10)+Tn(j/10)),s=-2..2,color=red),j=1..100):
S2:=seq(spacecurve(evalm(H(j/10)+Nn(j/10)),s=-2..2,color=red),j=1..100):
S3:=seq(spacecurve(evalm(H(j/10)+Bn(j/10)),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,t=0..10,color=blue):
display(d1, d2,d3,d4,axes=boxed,scaling=constrained);
|
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(plots): with(linalg):
|
Warning, new definition for norm
Warning, new definition for trace
> |
h:=[t*cos(t),2*t*sin(t),t]:
|
> |
`dh/dt`:=map(diff,h,t):
|
> |
`ds/dt`:=sqrt(dotprod(`dh/dt`,`dh/dt`)):
|
> |
`dh/ds`:=evalm(eval(`dh/dt`)/eval(`ds/dt`)):
|
> |
dT:=simplify(map(diff,eval(`dh/ds`),t)):
|
> |
n:=simplify(evalm(dT/eval(`ds/dt`))):
|
> |
rho:=j->evalf(1/(subs(t=j,simplify(sqrt(innerprod(n,n)))))):
|
> |
H:=i->subs(t=i,eval(h)):
|
> |
T:=p->subs(t=p,eval(`dh/ds`)):
|
> |
N:=q->subs(t=q,eval(n)/sqrt(innerprod(n,n))):
|
> |
Circ:=k->evalf(evalm(H(k/10)+rho(k/10)*(N(k/10)+(cos(theta)*T(k/10)+sin(theta)*N(k/10))))):
|
> |
S:=seq(spacecurve(Circ(k),theta=0..2*Pi),k=0..100):
|
> |
P1:=display(S,insequence=true,color=blue):
|
> |
P2:=spacecurve(h,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):
|
> |
h:=[t*cos(t),2*t*sin(t),t]:
|
> |
c1:=spacecurve(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(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.
To project onto the x-z plane we enter:
> |
T2:=transform((x,y,z)->[x,0,z]):
|
Practice
1. For the following curves, find and plot:
a) The unit tangent.
b) The unit normal
c) The unit binormal at the point iondicated.
i)
t=3.0
ii)
t=1
iii)
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.