Calculus IV with Maple
Copyright 2002, Dr. Jack Wagner
j.wagner@intelligentsearch.com

Lesson 6: Curvature

Topics:   Normal curvature, Gaussian curvature, second fundamental form, Dupin Indicatrix ,
Maple commands introduced : PrincipalNormal, BilinearForm, Curvature

Normal Curvature

Consider a surface,  in which we embed the image of the curve .

>	restart: with(VectorCalculus): with(LinearAlgebra):with(plots):

>	F := <u, v, u^2-v^2>;

>	plot3d(F, u=-2..2, v=-2..2, axes=normal, color=aquamarine, style=wireframe);

>	P := <t, cos(t), 0>;

>	C := eval(F, {u=t, v=cos(t)});

>	P0 := spacecurve(evalm(P), t=-2..2, color=blue, thickness=2):

>	P1 := plot3d(F, u=-2..2, v=-2..2, axes=normal, color=aquamarine, style=wireframe):

>	P2 := spacecurve(evalm(C), t=-2..2, color=red, thickness=2):

>	display(P0, P1, P2);

Now define a map,  . This is the unit normal to the surface and is a map to the unit sphere. Restricted to C it is a curve embedded in the surface of the unit sphere.

>	dF[u] := diff(F, u);

>	dF[v] := diff(F, v);

>	nn :=  (dF[u] &x dF[v]);

>	N := Normalize(nn, 2, inplace):

>	N := simplify(N) assuming real;

>	Nt := eval(N, {u=t, v=cos(t)});

The normal vector N(u,v) shown as a parametrized surface over u and v.

>	P3 := plot3d(N, u=-2..2, v=-2..2, color=aquamarine, axes=normal, style=wireframe, scaling=constrained, numpoints=2000):

>	P4 := spacecurve(evalm(Nt), t=-2..2, axes=normal, color=red, thickness=2, style=wireframe):

>	display(P3, P4);

>	dC := diff(C, t);

>	dC0 := eval(dC, t=0); #The tangent to C at t=0

>	C0 := eval(C, t=0);

>	P5 := spacecurve(evalm(C0+s*dC0), s=-1..1, color=blue, thickness=3):

>	display(P1, P2, P5);

The differential of the map N, dN, is a map from the tangent space of the F(u, v) surface to the tangent space of the unit sphere.  Since N.N=1 -> N'.N = 0, dN lives in the tangent space. dN =  +  

>	JN := Jacobian(N, [u, v, w]):

Let dNt be the map dN restricted to the curve C(t) = F(P(t)). dNt = N(C(t))' = N'(C(t))C'(t). This is the tangent to the image of the curve C(t) on the unit sphere.

>	JNt := simplify( eval(JN, {u=t, v=cos(t)}) );

>	dP := diff(P, t);

>	JNt0 := eval(JNt, t=0);

>	Tt := simplify(JNt.dP);

Evaluating at t=0; The image of the tangent to C at t=0.

>	Tt0 := eval(Tt, t=0);

At the point N(0, 1) corresponding to t=0, the image of dC0 under the map dNt  is the tangent to the image of C under the map N.

>	N0 := eval(Nt, t=0);

>	P6 := spacecurve(evalm(N0+s*Tt0), s=-1/2..1/2, color=blue, thickness=2):

>	display(P3, P4, P6);

Now consider the principal normal vector  to the embedded curve C.

PrincipalNormal

>	n := simplify( PrincipalNormal(C, t) );

>	n0 := eval(n, t=0);

Note that the principal normal, n, is not   a unit vector. n =   n , where n  is the unit normal.

>	P7 := spacecurve({evalm(C0+s*Tt0), evalm(C0+3*s*N0), evalm(C0+s*n0)}, s=-1..1, color=blue, thickness=3):

>	display(P1, P2, P7);

Define the normal curvature as the projection of n onto N.

>	kappa['n'] := N0.n0;

Because N is a unit vector, <N, n> = <N, n > =  cos(N, n ), where  is the usual curvature.  

Curvature

>	k := eval(Curvature(C, t), t=0);

>	cos(theta) := N0.Normalize(n0, 2);

>	simplify(k*cos(theta));

Second Fundamental Form

We now define a bilinear form   . Set dN = JNt. Let x (s) =    be a curve in the surface F(u, v), parametrized by arc length. Let (s)  be the unit tangent at s. Define:   ( (s)) =   -   <dN ( (s)), (s)   >      ( (s)) is called the second fundamental form of F at p. We will work at the point p  = x (0). Let N( x (s)) be the restriction of N to x (s).

Since N( x (0)) is normal to the tangent   x' (0) we have <N( x (0)), x' (0)>=0.

<N'( x (0). x' (0), x' (0)> + <N, x'' (0)> = 0 -> <N'( x (0)). x' (0), x' (0)> = - <N, x'' (0)> =

Implementing this in Maple:

>	II[p] := JNt0;

>	u := t: v := cos(t):

>	du := diff(u, t); dv := diff(v, t);

>	dx := <du, dv, 0>;

Note that at t=0 this is a unit  tangent vector.   

BilinearForm

>	B := -BilinearForm(dx, dx, JNt0);

>	Bt0 := eval(B, t=0);

Expand ( ) =   -   <dN ( ),   > = -   (   ).(  + ) = -   -     -   -   .   

Since N  is orthogonal to the tangent plane,    = 0 ->  = -  ->

( ) =  +2  +  

Set L=  , M =  , N =   The second fundamental form may then be expressed as:

( ) = L  + 2M  + N  

>	restart: with(VectorCalculus): with(LinearAlgebra):with(plots):

>	F := <u, v, u^2-v^2>;

>	dF[u] := diff(F, u); dF[v] := diff(F, v);

>	n := dF[u] &x dF[v];

>	NN := Normalize(n, 2, inplace):

>	NN := simplify(NN) assuming real;

>	dF[uu] := diff(dF[u], u);

>	dF[uv] := diff(dF[u], v);

>	dF[vv] := diff(dF[v], v);

>	L := NN.dF[uu];

>	M := NN.dF[uv];

>	N := NN.dF[vv];

>	u := t: v := cos(t):

>	du := diff(u, t); dv := diff(v, t);

>	II[p] := L*(du)^2+M*du*dv+N*(dv)^2;

>	eval(II[p], t=0);

>	restart:with(plots): with(VectorCalculus): with(LinearAlgebra):

The maximum and minimum curvatures are the eigenvalues of dN.  

>	F := <u, v, u^2-v^2>;

>	dF[u] := diff(F, u); dF[v] := diff(F, v);

>	n := dF[u] &x dF[v];

>	N := Normalize(n, 2, inplace):

>	N := simplify(N) assuming real;

>	dN := Jacobian(N, [u, v, w]):

>	dNt := subs({u=t, v=cos(t)}, dN):

>	dNt0 := subs(t=0, dNt);

>	E := Eigenvalues(dNt0);

>	kappa[1][t=0] := E[2];

>	kappa[2][t=0] := E[3];

      The eigenvectors of dN  are the principal directions.

>	E := Eigenvectors(dNt0, output='vectors');

>	e[1] := <0, -1/2, 1>;

>	e[2] := <1, 0, 0>;

>	C := subs({u=t, v=cos(t)}, F);

>	C0 := subs(t=0, C);

>	P1 := plot3d(F, u=-2..2, v=-2..2, color=aquamarine, style=wireframe, axes=normal):

>	P2 := spacecurve({evalm(C0+s*e[1]), evalm(C0+s*e[2])}, s=-1..1, color=red, thickness=2):

>	P3 := spacecurve(evalm(C), t=-2..2, color=blue, thickness=2):

>	display(P1, P2, P3);

Gaussian Curvature  

>	restart: with(LinearAlgebra): with(VectorCalculus):with(linalg):

The Gaussian curvature K  =  .  The mean curvature H =  .

In the coordinate basis of the ambient space , [u, v, w], let .

>	N := [n(u, v)[1], n(u, v)[2], n(u, v)[3]];

>	dN := Jacobian(N, [u, v, w]);

>	E := eigenvalues(dN);

>	K := expand(E[2]*E[3]);

>	H := (E[2]+E[3])/2;

Both K  and H  are independent of the z component of N, and are equal, respectively, to the determinant and trace of the minor of    

In the local  coordinate basis,  :

>	A := Matrix([[a[11], a[12]], [a[21], a[22]]]);

>	dN := A*Matrix([[x[u]], [x[v]]])=A.Matrix([[x[u]], [x[v]]]);

In this basis dN  is given by the A matrix.

The eigenvalues of the A matrix are:

>	e := Eigenvalues(A);

>	K := expand(e[1]*e[2]);

K is the determinant of the A matrix.

>	H := (e[1]+e[2])/2;

H is 1/2 Tr(A).

Since N  is a map from the surface, F ( u, v ), to the unit sphere, dN  is a map from the tangent space of   F ( u, v ) to the tangent space of the sphere. More specifically, let x  = [u(t), v(t)] be a curve embedded in the surface and  be a unit tangent vector at  point p .

( ) = < , >       

 ( ) = -<  dN ( ),  >

In matrix form, recollecting the designation of the coefficients of  and .

>	N := 'N':F := 'F':E := 'E':

>	II_[p] := <<L|M>, <M|N>>;

>	I_[p] := <<E|F>, <F|G>>;

>

A := A;

>	Eq1 := II_[p]=-A*I_[p];

>	invI_[p] := MatrixInverse(I_[p]);

>	Eq2 := A=simplify(-II_[p].invI_[p]);

This provides an explicit formula for the A matrix components.

Dupin Indicatrix

For unit vectors, v, in the basis formed by the orthogonal (but not orthonormal) curvature directions, the locus of points:

+/-1 =  is known as the Dupin indicatrix.

>	restart:with(plots): with(VectorCalculus): with(LinearAlgebra):

>	assume(u, real): assume(v, real): interface(showassumed=0):

>	F := <u, v, u^2-v^2>;

>	dF[u] := diff(F, u); dF[v] := diff(F, v);

>	n := dF[u] &x dF[v];

>	N := Normalize(n, 2, inplace);

>	N := simplify(N) assuming real;

>	dN := Jacobian(N, [u, v, w]):

>	dNt := subs({u=t, v=cos(t)}, dN):

>	dNt0 := subs(t=0, dNt);

>	E := Eigenvectors(dNt0, output='vectors');

>	e[1] := convert(Column(E, 1), Vector[row]);

>	e[2] := convert(Column(E, 2), Vector[row]);

>	V := VectorScalarMultiply(e[1], u)+VectorScalarMultiply(e[2], v);

>	eq1 := 1=eval(BilinearForm(V, V, dNt0));

>	eq2 := -1=eval(BilinearForm(V, V, dNt0));

Actually, the Dupin indicatrix is defined on normalized vectors but this would only complicate the arithmetic and add nothing to the picture.

>	implicitplot({eq1, eq2}, u=-5..5, v=-5..5, numpoints=1000);

This gives substance to the designation "hyperbolic point"  for those points on a surface with negative Gausian curvature.

>	Eg := Eigenvalues(dNt0);

>	K := Eg[2]*Eg[3];

In the case of a point with positive Gaussian curvature:

>	restart: with(LinearAlgebra): with(VectorCalculus): with(plots):

>	assume(u, real): assume(v, real):

>	interface(showassumed=0):

>	G := <u, v, u^2+v^2>;

>	dG[u] := diff(G, u); dG[v] := diff(G, v);

>	n := dG[u] &x dG[v];

>	N := Normalize(n, 2, inplace):

>	N := simplify(N) assuming real;

>	dN := Jacobian(N, [u, v, w]):

>	dNt := subs({u=t, v=cos(t)}, dN):

>	dNt0 := subs(t=0, dNt);

>	E := Eigenvectors(dNt0, output='vectors');

>	e[2] := <0, 1/2, 1>;

>	e[1] := <1, 0, 0>;

>	V := VectorScalarMultiply(e[1], u)+VectorScalarMultiply(e[2], v);

>	eq1 := 1=eval(BilinearForm(V, V, dNt0));

>	eq2 := -1=eval(BilinearForm(V, V, dNt0));

>	implicitplot({eq1, eq2}, u=-5..5, v=-5..5, numpoints=2000, scaling=constrained);

The designation "elliptic point"  for those points on a surface with positive Gaussian curvature seems appropriate.

>	Eg := Eigenvalues(dNt0);

>	K := Eg[2]*Eg[3];

Now we visualize the surface G(u, v), together with the Dupin indicatrix and the curvature directions at the point u=0, v=1.

First we recollect the curvature directions.

>	e[1] := e[1];

>	e[2] := e[2];

      Construct a third orthogonal axis.

>	e[3] := e[2] &x e[1];

     The matrix of basis transformation:

>	with(linalg): T := stackmatrix(e[1], e[2], e[3]);

Warning, the previous bindings of the names GramSchmidt and Wronskian have been removed and they now have an assigned value

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

>	whattype(T);

>	T := convert(T, Matrix);

The relationship between our local coordinate system and the ambient system is:  

>	eq := <<e_[1], e_[2], e_[3]>> = T*<<e[x], e[y], e[z]>>;

The matrix of coordinate transformation:

>	TTinv := Transpose(MatrixInverse(T));

Solve the equation of the indicatrix for one variable.

>	V := solve(eq2, v);

Define the two halves of the indicatrix as vectors.  

>	C1 := <u, V[1], 0>;

>	C2 := <u, V[2], 0>;

Transform the vector components with the proper transformation matrix.

>	Ct1 := TTinv.C1;

>	Ct2 := TTinv.C2;

Identify the point on the surface corresponding to u=0, v=1.

>	G0 := subs({u=0, v=1}, G);

Plot the indicatrix and the coordinate directions at the indicated pointon the surface, together with a plot of the surface.

>	P1 := spacecurve({evalm(G0+Ct1), evalm(G0+Ct2)}, u=-5..5, numpoints=1000, color=blue, axes=normal, scaling=constrained):

>	P2 := spacecurve({evalm(G0+s*e[1]), evalm(G0+s*e[2])}, s=-2..2, color=red, thickness=2):

>	P3 := plot3d(G, u=-1..1, v=0..1.5, color=aquamarine, style=wireframe):

>	display(P1, P2, P3);