Geometry of Curves
by John Oprea
This worksheet contains some basic applications of MAPLE to the differential geometry of curves. A general reference is J. Oprea, Differential Geometry and its Applications .
In particular, there are procedures for computing the curvature and torsion of curves, and for determining a curve solely from its curvature and/or torsion.
The first set of procedures are for dot product, the length (or norm) of a vector and for the cross product.
> dp := proc(X,Y)
> X[1]*Y[1]+X[2]*Y[2]+X[3]*Y[3];
> end:
> nrm := proc(X)
> sqrt(dp(X,X));
> xp := proc(X,Y)
> local a,b,c;
> a := X[2]*Y[3]-X[3]*Y[2];
> b := X[3]*Y[1]-X[1]*Y[3];
> c := X[1]*Y[2]-X[2]*Y[1];
> [a,b,c];
These procedures compute curvature and torsion. These functions assume that the curve is given parametrically in three coordinates. If the curve is a plane curve in the xy-plane, say, then make the third coordinate zero.
> curv:=proc(alpha)
> local alphap,alphapp;
> alphap:=diff(alpha,t);
> alphapp:=diff(alphap,t);
> RETURN(kappa=simplify(nrm(xp(alphap,alphapp))/nrm(alphap)^3,symbolic));
> tor:= proc(alpha)
> local alphap,alphapp,alphappp;
> alphappp:=diff(alphapp,t);
> RETURN(tau=simplify(dp(xp(alphap,alphapp),alphappp)/nrm(xp(alphap,alphapp))^2,symbolic));
Here's an example. Take a helix.
> hel:=[a*cos(t),a*sin(t),b*t];
> curv(hel);
> tor(hel);
> hel1:=subs({a=4,b=2},hel);
> curv(hel1);
> tor(hel1);
Before special types of plots can be made, the following command must be issued.
> with(plots):
The following plot command graphs the helix in 3-space and colors it red.
> spacecurve(hel1,t=0..10*Pi,scaling=constrained,thickness=3,color=red);
To plot plane curves, we can use the plot command. For example, take the
witch of Agnesi.
> witch:=[2*tan(t),2*cos(t)^2,0];
> curv(witch);
You cannot say plot(witch ...) since you've used three coordinates here.
You must say (in THIS form)
> plot([2*tan(t),2*cos(t)^2,t=-1.3..1.3],scaling=constrained,axes=framed,color=red);
or use the first two coordinates of the witch as follows
> plot([witch[1],witch[2],t=-1.3..1.3],scaling=constrained,axes=framed,color=red);
Notice that the square brackets must go around the t=... as well as the parametrization. Try using options like `constrained` or `color` or `axes` or `thickness`. It is a basic theorem that a plane curve may be recreated from its curvature alone. While the formula involves integrals which can only rarely be solved, we can still plot the resulting curve. Given curvature , the formula is
> beta(s)=[Int(cos(theta(u)),u=0..s),Int(sin(theta(u)),u=0..s)];
where
> theta(u)=Int(kappa(t),t=0..u);
By the Fundamental Theorem of Calculus, we can transform these integrals into a system of differential equations to be solved numerically and plotted. The result will be the unit speed curve with the specified curvature. First, let's take itself.
> sys:=diff(theta(s),s)=s, diff(b1(s),s)=cos(theta(s)), diff(b2(s),s)=sin(theta(s));
> p:=dsolve({sys,theta(0)=0,b1(0)=0,b2(0)=0},{theta(s),b1(s),b2(s)},type=numeric):
> test:=odeplot(p,[b1(s),b2(s)],-5..5,numpoints=200):
> display(test,view=[-2..2,-2..2]);
As a second example, take curvature=1/(1+s^2).
> sys:=diff(theta(s),s)=1/(1+s^2), diff(b1(s),s)=cos(theta(s)), diff(b2(s),s)=sin(theta(s));
> odeplot(p,[b1(s),b2(s)],-8..8,numpoints=200,scaling=constrained);
What is this curve? Can you prove it?
Now let's write a procedure to take an input curvature and output the curve. Other inputs give a bound on the curves parameter(a,b) and bounds on the viewing region (c,d,f,g).
> recreate:=proc(kap,a,b,c,d,f,g)
> local sys,b1,b2,p,theta,pl;
> sys:=diff(theta(s),s)=kap(s), diff(b1(s),s)=cos(theta(s)), diff(b2(s),s)=sin(theta(s));
> pl:=odeplot(p,[b1(s),b2(s)],a..b,numpoints=400,thickness=1,axes=framed, color=red):
> display(pl,view=[c..d,f..g],scaling=constrained);
To specify a curvature function, give the parameter (t say) and use the hyphen and greater than sign to construct an arrow -> to the effect on t. Try these.
> kap1:=t->t;
> kap2:=t->t^2;
> kap3:=t->exp(t);
> kap4:=t->sin(t);
> kap5:=t->sin(t)*t;
> kap6:=t->sin(t)*t^2;
> kap7:=t->1/(1+t^2);
> recreate(kap5,-8,8,-2,2,0,3);
> recreate(kap7,-10,10,-4,4,0,10);
We can also recreate curves in space from their curvature and torsion via the Frenet (differential) equations. The following procedure solves the nine Frenet equations together with the three equations defining the unit tangent vector of the curve and allows you to zoom in by defining a viewing region. The inputs c,d,e,f,g,h give ranges for x,y and z respectively.
> recreate3dview:=proc(kap,ta,a,b,c,d,e,f,g,h)
> local sys,p,alph1,alph2,alph3,T1,T2,T3,N1,N2,N3,B1,B2,B3,pl;
> sys:=diff(alph1(s),s)=T1(s),diff(alph2(s),s)=T2(s),diff(alph3(s),s)=T3(s),
> diff(T1(s),s)=kap(s)*N1(s), diff(T2(s),s)=kap(s)*N2(s),
> diff(T3(s),s)=kap(s)*N3(s),diff(N1(s),s)=-kap(s)*T1(s)+ta(s)*B1(s),
> diff(N2(s),s)=-kap(s)*T2(s)+ta(s)*B2(s),diff(N3(s),s)=-kap(s)*T3(s)+ta(s)*B3(s),
> diff(B1(s),s)=-ta(s)*N1(s), diff(B2(s),s)=-ta(s)*N2(s),
> diff(B3(s),s)=-ta(s)*N3(s); print(sys);
> p:=dsolve({sys, alph1(0)=0,alph2(0)=0,alph3(0)=0,T1(0)=1,T2(0)=0,T3(0)=0,
> N1(0)=0,N2(0)=1,N3(0)=0,B1(0)=0,B2(0)=0,B3(0)=1},
> {alph1(s),alph2(s),alph3(s),T1(s),T2(s),T3(s),N1(s),N2(s),N3(s),
> B1(s),B2(s),B3(s)},type=numeric):
> pl:=odeplot(p,[alph1(s),alph2(s),alph3(s)],a..b,numpoints=200,thickness=1,
> axes=framed, color=red):
> display(pl,scaling=constrained,view=[c..d,e..f,g..h]);
Here are some examples. What is the first curve?
> recreate3dview(1/5,1/10,0,200,0,25,0,20,0,40);
> kap3d1:=t->t;
> recreate3dview(kap3d1,1/10,0,10,0,2,0,2,0,0.5);
> tau3d1:=t->t/10;
> recreate3dview(kap3d1,tau3d1,0,10,0,2,0,2,0,1);
Here is the general formula for a hypotrochoid:
> hypogen:=[(a-b)*cos(t)+h*cos((a-b)*t/b),(a-b)*sin(t)-h*sin((a-b)*t/b),0];
> curv(subs({a=1,b=1/3,h=2/3},hypogen));
As a test of our procedure recreate, let's use the curvature formula to recreate a hypotrochoid and the Witch and then compare the results to ones drawn directly.
> hypocurv:=t->-3/2*(-7-6*cos(t)+8*cos(t)^3)/((-5+16*cos(t)^3-12*cos(t))*sqrt(5-16*cos(t)^3+12*cos(t))):
> plot([subs({a=1,b=1/3,h=2/3},hypogen)[1],subs({a=1,b=1/3,h=2/3},hypogen)[2],t=0..2*Pi]);
> recreate(hypocurv,0,2*Pi,-1,1,-1.5,0.2);
> witchcurv:=t->-(4*cos(t)^2-3)*cos(t)^4/((-4*cos(t)^6-1+4*cos(t)^8)*sqrt(4*cos(t)^6+ 1-4*cos(t)^8));
> recreate(witchcurv,-3,3,-1,1,0,0.4);
