• Define and Plot the Curve
• Velocity, Acceleration and Jerk
• Speed and Arclength
• Frenet Frame
• Plot the Frame
• Curvature and Torsion
• Tangential and Normal Acceleration

Frenet Analysis of Curves

Using Maple and the vec_calc Package

In this worksheet we will see how Maple and the vec_calc package can be used to analyse a parametrized curve. As an example, we will consider:

* The Winding Line on a Torus

which may be parametrized by

Here is the big radius, is the small radius and and are the angular rates.

The quantities we will compute are:

* the velocity

* the acceleration

* the jerk

* the speed

* the arclength

* the unit tangent vector

* the unit binormal vector

* the unit principal normal vector

* the curvature =

* the torsion =

* the tangential acceleration =

* the normal acceleration =

In addition we will animate the drawing of the curve .

We start by initializing the vec_calc package:

> restart;

> libname:=libname,"C:/mylib/vec_calc7":

> with(vec_calc): vc_aliases:

> with(linalg):with(student):with(plots):

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the name changecoords has been redefined

>

Define and Plot the Curve

We enter the function using the makefunction command

> r:=MF(t, [(p+q*cos(n*t))*cos(m*t), (p+q*cos(n*t))*sin(m*t), q*sin(n*t)]);

We choose values for the parameters: As goes from to , tells how many times the curve wraps around the circle of radius and tells how many times the curve wraps around the circle of radius . If the ratio is rational, the curve closes. If is irrational, the curve does not close and has infinite length.

> p:=1; q:=1/2; m:=1; n:=8;

We plot the curve: (Notice that the op command is needed to strip off the square brackets from the result of .)

> spacecurve([op(r(t)),t=0..2*Pi], numpoints=193, scaling=constrained, orientation=[0,45], color=blue);

Try the following:

* Click in the plot and drag the mouse. This will give you a better idea about the 3-dimensional shape of the curve.

* Change the numbers and . In particular, try , , or try , , or try , . For the last two, the curve does not close.

We can also animate the curve being drawn:

> frames:=24;

> display([seq(spacecurve([op(r(t)), t=0..i*2*Pi/frames], numpoints=8*i, color=blue), i=1..frames)], orientation=[0,45], scaling=constrained, insequence=true);

Click in the plot. Then click on the loop and play buttons. Also try dragging the mouse while the plot is running.

We will add to these plots later.

Before going on, we remove the values of the parameters so the remaining quantites are computed in general whenever possible.

> p:='p'; q:='q'; m:='m'; n:='n';

We will reassign these parameters whenever necessary.

>

Velocity, Acceleration and Jerk

The velocity is the derivative of the position: (Notice how we can handle the curve as either a vector of expressions or as a vector of functions.)

> diff(r(t),t);

> v:=D(r);

> v(t);

The acceleration is the derivative of the velocity:

> a:=D(v);

> a(t);

The jerk is the derivative of the acceleration:

> j:=D(a);

> j(t);

There are also vec_calc commands to directly compute the velocity, acceleration and jerk:

> v2:=Cv(r);

> a2:=Ca(r);

> j2:=Cj(r);

We check the answers are the same:

> v(t)-v2(t);

> a(t)-a2(t);

> j(t)-j2(t);

>

Speed and Arclength

The speed is the length of the velocity:

> speed:=simplify(len(v(t)));

The arclength is the integral of the speed: ( Maple is unable to compute the exact value using the value command.)

> Int(speed,t=0..2*Pi); #arclength:=value(%);

So we give values to the parameters and use evalf to get a decimal approximation:

> p:=1; q:=1/2; m:=1; n:=8;

> Int(speed,t=0..2*Pi); arclength:=evalf(%);

Notice that since , , the arclength is a little more than times the circumference of a circle of radius :

> 8*2*Pi*q; evalf(%);

There is also a vec_calc command to directly compute the arclength integral as a function of its two endpoints:

> L:=CL(r);

> L(0,2*Pi); evalf(%);

Once again, we unassign the parameters:

> p:='p'; q:='q'; m:='m'; n:='n';

>

Frenet Frame

The vectors , and form a right handed triplet of unit vectors called the Frenet frame for the curve.

The unit tangent vector is the unit vector in the direction of the velocity:

> T:=v(t)/speed;

The unit binormal is the unit vector perpendicular to the instantaneous plane of the motion in the direction of x . We first need

> vxa:=simplify(v(t) &x a(t));

> lvxa:=simplify(len(vxa));

Then the unit binormal is

> B:=vxa/lvxa;

Since , and form a right handed triplet, in that order, we have that the unit principal normal is: (The evall command multiplies the scalars into the vectors.)

> N:=simplify(evall(B) &x evall(T));

There are also vec_calc commands to directly compute , and as vector functions:

> T2:=CT(r);

> B2:=CB(r);

> N2:=CN(r);

We again check the answers are the same:

> simplify(evall(T-T2(t)));

> evall(B-B2(t));

Maple is unable to show the difference of the vectors is zero. So it needs some help. We show the square of each component is the same by showing the ratio is one:

> N[1]^2/N2(t)[1]^2;

> N[2]^2/N2(t)[2]^2;

> N[3]^2/N2(t)[3]^2;

>

Plot the Frame

We want to add the Frenet frame to our plot and animation of the curve.

We first redefine the parameters.

> p:=1; q:=1/2; m:=1; n:=8;

Next we replot the curve but stop at a random endpoint:

> tfinal:=1.2*Pi;

> crvplot:=spacecurve([op(r(t)),t=0..tfinal], numpoints=193, scaling=constrained, orientation=[0,45], color=blue): crvplot;

Now we plot the tangent vector using a dot-to-dot plot, but locate the vector at the tip of the curve: (Note: We need to use evall to add the two vectors.)

> Tplot:=spacecurve([r(tfinal),evall(r(tfinal)+T2(tfinal))], color=red): Tplot;

In addition, we plot the normal and the binormal and display all three with the original curve:

> Nplot:=spacecurve([r(tfinal),evall(r(tfinal)+N2(tfinal))], color=cyan):

> Bplot:=spacecurve([r(tfinal),evall(r(tfinal)+B2(tfinal))], color=magenta):

> display({crvplot,Tplot,Nplot,Bplot});

We are finally ready to animate the curve:

> frames:=50;

> plotlist:=NULL:

> for i from 1 to frames do
tfinal:=i*2*Pi/frames;
crvplot:=spacecurve([op(r(t)),t=0..tfinal], numpoints=193, scaling=constrained, orientation=[0,45], color=blue):
Tplot:=spacecurve([r(tfinal),evall(r(tfinal)+T2(tfinal))], color=red):
Nplot:=spacecurve([r(tfinal),evall(r(tfinal)+N2(tfinal))], color=cyan):
Bplot:=spacecurve([r(tfinal),evall(r(tfinal)+B2(tfinal))], color=magenta):
plotlist:=plotlist,display({crvplot,Tplot,Nplot,Bplot});
od:

> display([plotlist], orientation=[0,45], scaling=constrained, insequence=true);

Before going on, we again remove the values of the parameters.

> p:='p'; q:='q'; m:='m'; n:='n';

>

Curvature and Torsion

The curvature measures the rate at which the curve is turning. There are a couple ways to compute the curvature. We do two and compare them:

> kappa:=(a(t) &. N)/speed^2;

> Kappa:=simplify(lvxa/speed^3);

> kappa-Kappa;

The torsion measures the rate at which the plane of the curve is turning.

> tau:=(vxa &. j(t))/lvxa^2;

There are also vec_calc commands to directly compute the curvature and torsion as functions

> kappa2:=Ck(r);

> tau2:=Ct(r);

We again check these agree: Again Maple had trouble simplifying the difference of the curvatures, so we compute the ratio of the squares:

> simplify(kappa^2/kappa2(t)^2);

> simplify(tau-tau2(t));

>

Tangential and Normal Acceleration

The acceleration can be decomposed into components parallel and perpendicular to the direction of motion: The coefficients are the tangential and normal accelerations. We compute each of them in two ways and compare.

> aT:=a(t) &. evall(T);

> AT:=simplify(diff(speed,t));

> aT-AT;

> aN:=a(t) &. evall(N);

> AN:=kappa*speed^2;

> aN-AN;

Again there vec_calc commands to directly compute the tangential and normal accelerations.

> aT2:=CaT(r);

> aN2:=CaN(r);

Again we check:

> aT-aT2(t);

> aN^2-aN2(t)^2;simplify(%);

>