crv.mws
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(%);
>