Math 276

Calculus III

Spring 2003

Dr. Constant J. Goutziers

Department of Mathematics, Computer Science and Statistics

goutzicj@oneonta.edu

Lesson 8

Arclength, Curvature

the Principal Normal and Binormal vectors

Initializations

>    restart;
with(linalg):
with(oneonta):
with(plots):

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

Warning, the name changecoords has been redefined

8.1  Arclength

The arclength of a space curve r = r(t) equals Int(abs(diff(r(t),t)),t = a .. b) .

Examples

Example 8.1.1
Compute the arclength of the curve r(t) = [cos(7*t), exp(-t), ln(10*t+1)] , 0 `` <= ``  t `` <= `` 5.

Define the curve and the integral representing the arclength.

>    r:=vector([cos(t), exp(-t), ln(10*t+1)]);

r := vector([cos(t), exp(-t), ln(10*t+1)])

>    rp:=map(diff, r, t);

rp := vector([-sin(t), -exp(-t), 10/(10*t+1)])

>    arclength:=Int(realnorm(rp), t=0..5);

arclength := Int((sin(t)^2+exp(-t)^2+100/(10*t+1)^2)^(1/2),t = 0 .. 5)

Of course it is not possible to give closed form representation of this integral.  In stead we produce a decimal approximation.

>    arclength_decimal:=evalf(arclength);

arclength_decimal := 5.930146803

8.2  Curvature

The curvature kappa  of a spacecurve is defined as the norm of   diff(T(s),s) .  Because T has constant length, we know that diff(T(t),t)  and therefore diff(T(s),s)  is perpendicular to T.  This indicates that kappa  truly measures the rate of change of the direction of the tangent vector to the curve.  We compute diff(T(s),s)  as   diff(T(t),t)/diff(s(t),t)   , where diff(s(t),t) = abs(diff(r(t),t)) .

Examples

Example 8.2.1
Compute the curvature of r(t) = [t, t^2, t^3] .

Just follow the theory described in the introduction of curvature.

>    r:=vector([t, t^2, t^3]);

r := vector([t, t^2, t^3])

>    rp:=map(diff, r, t);

rp := vector([1, 2*t, 3*t^2])

>    ds_dt:=realnorm(rp);

ds_dt := (1+4*t^2+9*t^4)^(1/2)

>    T:=evalm(rp/ds_dt);

T := vector([1/((1+4*t^2+9*t^4)^(1/2)), 2/(1+4*t^2+9*t^4)^(1/2)*t, 3/(1+4*t^2+9*t^4)^(1/2)*t^2])

>    dT_ds:=map(simplify, evalm(map(diff, T, t)/ds_dt));

dT_ds := vector([-2*t*(2+9*t^2)/(1+4*t^2+9*t^4)^2, -2*(9*t^4-1)/(1+4*t^2+9*t^4)^2, 6*t*(2*t^2+1)/(1+4*t^2+9*t^4)^2])

>    kappa:=simplify(realnorm(dT_ds));

kappa := 2*((9*t^4+9*t^2+1)/(1+4*t^2+9*t^4)^3)^(1/2)

>    kappa:=simplify(kappa, symbolic);

kappa := 2/(1+4*t^2+9*t^4)^(3/2)*(9*t^4+9*t^2+1)^(1/2)

>   

Example 8.2.2
Find a formula for the curvature of r(t) = [t*ln(t), sin(3*t)] .   Sketch the curve and its curvature in one picture.

>    r:=vector([t*ln(t), sin(3*t)]);

r := vector([t*ln(t), sin(3*t)])

>    rp:=map(diff, r, t);

rp := vector([ln(t)+1, 3*cos(3*t)])

>    ds_dt:=realnorm(rp);

ds_dt := ((ln(t)+1)^2+9*cos(3*t)^2)^(1/2)

>    T:=map(simplify, evalm(rp/ds_dt));

T := vector([1/(ln(t)^2+2*ln(t)+1+9*cos(3*t)^2)^(1/2)*(ln(t)+1), 3/(ln(t)^2+2*ln(t)+1+9*cos(3*t)^2)^(1/2)*cos(3*t)])

>    dT_ds:=map(simplify, evalm(map(diff, T, t)/ds_dt));

dT_ds := vector([9*cos(3*t)*(3*ln(t)*sin(3*t)*t+3*sin(3*t)*t+cos(3*t))/t/(ln(t)^4+4*ln(t)^3+6*ln(t)^2+18*ln(t)^2*cos(3*t)^2+4*ln(t)+36*ln(t)*cos(3*t)^2+1+18*cos(3*t)^2+81*cos(3*t)^4), -3*(cos(3*t)*ln(t...
dT_ds := vector([9*cos(3*t)*(3*ln(t)*sin(3*t)*t+3*sin(3*t)*t+cos(3*t))/t/(ln(t)^4+4*ln(t)^3+6*ln(t)^2+18*ln(t)^2*cos(3*t)^2+4*ln(t)+36*ln(t)*cos(3*t)^2+1+18*cos(3*t)^2+81*cos(3*t)^4), -3*(cos(3*t)*ln(t...

>    kappa:=simplify(realnorm(dT_ds));

kappa := 3*(-(-9*t^2*ln(t)^2+9*t^2*ln(t)^2*cos(3*t)^2-18*t^2*ln(t)+18*t^2*ln(t)*cos(3*t)^2-6*ln(t)*cos(3*t)*sin(3*t)*t-9*t^2+9*t^2*cos(3*t)^2-6*cos(3*t)*sin(3*t)*t-cos(3*t)^2)/t^2/(ln(t)^2+2*ln(t)+1+9*...

Observe that in order to synchronize the display of the curvature and the curve itself, we need to use equal x-scaling in both representations.

>    curvature:=plot([r[1], kappa, t=1..5]):

>    curve:=plot([r[1], r[2], t=1..5], color=blue):

>    with(plots):

>    display([curvature, curve]);

[Maple Plot]

>   

Example 8.2.3
An alternative formula for the curvature of a space curve in R^3  is: kappa = abs(crossprod(diff(r(t),t),diff(r(t),`$`(t,2))))/(abs(diff(r(t),t))^3) .
Use this formula on the space curve given in example 8.2.1.

>    r:=vector([t, t^2, t^3]);

r := vector([t, t^2, t^3])

>    kappa:=realnorm(crossprod(map(diff, r, t), map(diff, r, t$2)))/realnorm(map(diff, r, t))^3;

kappa := 2/(1+4*t^2+9*t^4)^(3/2)*(9*t^4+9*t^2+1)^(1/2)

>   

Example 8.2.4
The formula in example 8.2.3 can readily be used to find a simple expression for the curvature of  y=f(x).

>    r:=vector([x, f(x), 0]);

r := vector([x, f(x), 0])

>    kappa:=realnorm(crossprod(map(diff, r, x), map(diff, r, x$2)))/realnorm(map(diff, r, x))^3;

kappa := (diff(f(x),`$`(x,2))^2)^(1/2)/(1+diff(f(x),x)^2)^(3/2)

When f is a real valued function the numerator of this expression can be replaced by abs(diff(f(x),`$`(x,2))) .

8.3  The principal Normal and the Binormal vector

The principal Normal vector N is defined as the unit vector in the direction of diff(T(t),t) .  The Binormal vector B is defined as the crossproduct of T and N.

Examples

Example 8.3.1
Let r(t) = [t, ln(cos(t))]  , -Pi/2  < t < Pi/2 .  Compute N and B.

>    r:=vector([t, ln(cos(t))]);

r := vector([t, ln(cos(t))])

>    rp:=map(diff, r, t);

rp := vector([1, -sin(t)/cos(t)])

>    T:=evalm(rp/realnorm(rp));

T := vector([1/((1+sin(t)^2/cos(t)^2)^(1/2)), -1/(1+sin(t)^2/cos(t)^2)^(1/2)*sin(t)/cos(t)])

Because cos(t) > 0 on the interval involved, this expression can be simplified under that assumption.

>    T:=map(simplify, T, symbolic);

T := vector([cos(t), -sin(t)])

>    N:=map(simplify, evalm(map(diff, T, t)/realnorm(map(diff, T, t))));

N := vector([-sin(t), -cos(t)])

In order to find B, we must first represent T and N as vectors in 3-space.

>    T:=[T[1], T[2], 0];
N:=[N[1], N[2], 0];

T := [cos(t), -sin(t), 0]

N := [-sin(t), -cos(t), 0]

>    B:=map(simplify, crossprod(T,N));

B := vector([0, 0, -1])

>   

Example 8.3.2
Find the equation of the osculating circle at the point P = [1, 1/2]  to the parabola y = x^2/2 .  Graph the osculating circle and the parabola in one picture.

First code the parabola as a space curve, then compute T, N and   kappa .

>    r:=[x, x^2/2,0];
rp:=diff(r, x);
rpp:=diff(r, x$2);

r := [x, 1/2*x^2, 0]

rp := [1, x, 0]

rpp := [0, 1, 0]

>    T:=evalm(rp/realnorm(rp));

T := vector([1/((1+x^2)^(1/2)), 1/(1+x^2)^(1/2)*x, 0])

>    Tp:=map(diff, T, x);
N:=map(simplify, evalm(Tp/realnorm(Tp)), symbolic);

Tp := vector([-1/(1+x^2)^(3/2)*x, -1/(1+x^2)^(3/2)*x^2+1/((1+x^2)^(1/2)), 0])

N := vector([-1/(1+x^2)^(1/2)*x, 1/((1+x^2)^(1/2)), 0])

>    kappa:=simplify(realnorm(crossprod(rp, rpp))/realnorm(rp)^3);

kappa := 1/((1+x^2)^(3/2))

Compute kappa  and N at the point P.

>    kappa_P:=subs(x=1, kappa);

kappa_P := 1/4*2^(1/2)

>    N_P:=subs(x=1, evalm(N));

N_P := vector([-1/2*2^(1/2), 1/2*2^(1/2), 0])

Compute the radius a  of the osculating circle.

>    a:=1/kappa_P;

a := 2*2^(1/2)

Compute the center C of the osculating circle.

>    C:=evalm(subs(x=1, r)+a*N_P);

C := vector([-1, 5/2, 0])

Compute the equation of the osculating circle.

>    eq:=(x-C[1])^2+(y-C[2])^2=a^2;

eq := (x+1)^2+(y-5/2)^2 = 8

Plot the osculating circle and the parabola in one picture using the implicit plot routine.

>    p1:=implicitplot(y=r[2], x=-10..10, y=-10..10, color=red, numpoints=4900):
p2:=implicitplot(eq, x=-10..10, y=-10..10, color=blue, numpoints=4900):
display([p1, p2], thickness=4, scaling=constrained);

[Maple Plot]

A higher quality picture can be obtained by parametrizing the osculating circle and using parametric plots instead.

>    c:=[-1+a*cos(t), 5/2+a*sin(t), 0];

c := [-1+2*2^(1/2)*cos(t), 5/2+2*2^(1/2)*sin(t), 0]

>    p1:=plot([c[1], c[2], t=0..2*Pi], color=blue):
p2:=plot(x^2/2, x=-4..4):
display([p1,p2], scaling=constrained);

[Maple Plot]

Observe that the osculating circle actually crosses the parabola at the point of tangency. This occurs because, for -1 `` < `` x < 1 , the curvature of the parabola is larger than the curvature of the osculating circle and for 1 < x ,  the curvature of the parabola is less than the curvature of the osculating circle.

It may be illustrative to show the ray from the center of the osculating circle to the point of tangency.  See the picture below.

>    P:=subs(x=1, r);
rr:=evalm(P+t*N_P);

P := [1, 1/2, 0]

rr := vector([1-1/2*t*2^(1/2), 1/2+1/2*t*2^(1/2), 0])

>    p3:=plot([rr[1], rr[2], t=0..a], color=magenta):

>    display([p1, p2, p3], scaling=constrained);

[Maple Plot]

>