Application Center - Maplesoft

App Preview:

Determining a parametrized curve from its curvature and torsion

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


 

Curves.mws

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));

> end:

> 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];

> end:

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));

> end:

> tor:= proc(alpha)

> local alphap,alphapp,alphappp;

> alphap:=diff(alpha,t);

> alphapp:=diff(alphap,t);

> alphappp:=diff(alphapp,t);

> RETURN(tau=simplify(dp(xp(alphap,alphapp),alphappp)/nrm(xp(alphap,alphapp))^2,symbolic));

> end:

Here's an example. Take a helix.

> hel:=[a*cos(t),a*sin(t),b*t];

hel := [a*cos(t), a*sin(t), b*t]

> curv(hel);

kappa = a/(b^2+a^2)

> tor(hel);

tau = b/(b^2+a^2)

> hel1:=subs({a=4,b=2},hel);

hel1 := [4*cos(t), 4*sin(t), 2*t]

> curv(hel1);

kappa = 1/5

> tor(hel1);

tau = 1/10

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);

[Maple Plot]

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];

witch := [2*tan(t), 2*cos(t)^2, 0]

> curv(witch);

kappa = -(4*cos(t)^2-3)*cos(t)^4/((-4*cos(t)^6-1+4*...

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);

[Maple Plot]

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);

[Maple Plot]

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 kappa , the formula is

> beta(s)=[Int(cos(theta(u)),u=0..s),Int(sin(theta(u)),u=0..s)];

beta(s) = [Int(cos(theta(u)),u = 0 .. s), Int(sin(t...

where

> theta(u)=Int(kappa(t),t=0..u);

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 kappa = s itself.

> sys:=diff(theta(s),s)=s, diff(b1(s),s)=cos(theta(s)), diff(b2(s),s)=sin(theta(s));

sys := diff(theta(s),s) = s, diff(b1(s),s) = cos(th...

> 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]);

[Maple Plot]

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));

sys := diff(theta(s),s) = 1/(1+s^2), diff(b1(s),s) ...

> p:=dsolve({sys,theta(0)=0,b1(0)=0,b2(0)=0},{theta(s),b1(s),b2(s)},type=numeric):

> odeplot(p,[b1(s),b2(s)],-8..8,numpoints=200,scaling=constrained);

[Maple Plot]

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));

> p:=dsolve({sys,theta(0)=0,b1(0)=0,b2(0)=0},{theta(s),b1(s),b2(s)},type=numeric):

> 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);

> end:

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;

kap1 := proc (t) options operator, arrow; t end pro...

> kap2:=t->t^2;

kap2 := proc (t) options operator, arrow; t^2 end p...

> kap3:=t->exp(t);

kap3 := exp

> kap4:=t->sin(t);

kap4 := sin

> kap5:=t->sin(t)*t;

kap5 := proc (t) options operator, arrow; sin(t)*t ...

> kap6:=t->sin(t)*t^2;

kap6 := proc (t) options operator, arrow; sin(t)*t^...

> kap7:=t->1/(1+t^2);

kap7 := proc (t) options operator, arrow; 1/(1+t^2)...

> recreate(kap5,-8,8,-2,2,0,3);

[Maple Plot]

> recreate(kap7,-10,10,-4,4,0,10);

[Maple Plot]

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]);

> end:

Here are some examples. What is the first curve?

> recreate3dview(1/5,1/10,0,200,0,25,0,20,0,40);

diff(alph1(s),s) = T1(s), diff(alph2(s),s) = T2(s),...
diff(alph1(s),s) = T1(s), diff(alph2(s),s) = T2(s),...

[Maple Plot]

> kap3d1:=t->t;

kap3d1 := proc (t) options operator, arrow; t end p...

> recreate3dview(kap3d1,1/10,0,10,0,2,0,2,0,0.5);

diff(alph1(s),s) = T1(s), diff(alph2(s),s) = T2(s),...
diff(alph1(s),s) = T1(s), diff(alph2(s),s) = T2(s),...

[Maple Plot]

> tau3d1:=t->t/10;

tau3d1 := proc (t) options operator, arrow; 1/10*t ...

> recreate3dview(kap3d1,tau3d1,0,10,0,2,0,2,0,1);

diff(alph1(s),s) = T1(s), diff(alph2(s),s) = T2(s),...
diff(alph1(s),s) = T1(s), diff(alph2(s),s) = T2(s),...

[Maple Plot]

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];

hypogen := [(a-b)*cos(t)+h*cos((a-b)*t/b), (a-b)*si...

> curv(subs({a=1,b=1/3,h=2/3},hypogen));

kappa = -3/2*(-7-6*cos(t)+8*cos(t)^3)/((-5+16*cos(t...

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]);

[Maple Plot]

> recreate(hypocurv,0,2*Pi,-1,1,-1.5,0.2);

[Maple Plot]

> 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));

witchcurv := proc (t) options operator, arrow; -(4*...

> recreate(witchcurv,-3,3,-1,1,0,0.4);

[Maple Plot]