Geodesic Curvature

Here is some Mathematica code to compute the geodesic curvature of a curve. It is the computation done in detail class last week. The example computes the geodesic curvature of a circle of lattiude curve in the unit sphere. Both the geodesic curvature as well as the normal curvature are signed like the curvature of a curve in the two dimensional plane. The actual curvature is not signed (meaning that it is non-negative). Think of the normal curvature as the scalar projection of the acceleration onto the normal vector and the geodesic curvature onto the crossed producte of the velocity with the normal vector. The formula is simpler if one has a perametrization with constant speed as then { x',n, x'x n } form an orthonormal frame along the curve similarly as the Frenet frame. But it is not the Frenet frame! The surface defines us what is "up" and not the curve itself. Also, in the geodesic case, geodesic curvature zero and the curvature is the absolute value of the normal curvature. For constant speed situations we have:

κg = (n x x'). x''/|x'|^3 κn = n. x''/|x'|^2 κ = |x''|/|x'|^2

Since the absolute values are lengths of a right angle triangle, Pythagoras gives:

κg2 + κn2 = κ2

In textbooks, you find these formulas in the case of parametrization with arc length where the |x'|=1. It is convenient to get the constant speed case as reparametrization can be a pain in the neck, even if it is only done by a scalar.

(* Example code to compute the geodesic curvature, O. Knill, Math 136 *) r={Sin[s] Cos[t],Sin[s] Sin[t],Cos[s]}; (* unit sphere *) ds=D[r,s]; dt=D[r,t]; ddt=D[r,{t,2}]; (* t -> r(s0,t) latitude*) n=Simplify[Cross[dt,ds]]; n = n/Sqrt[n.n]; (* normal to surface *) speed=Sqrt[dt.dt]; (* Speed *) normal=Simplify[ddt.n/speed^2] (* Normal curvature *) geodesic=Simplify[Cross[n,dt].ddt/speed^3] (* Geodesic curvature *) curvature=Simplify[Sqrt[ddt.ddt]]/speed^2 (* usual curvature *) Simplify[curvature^2 == normal^2+geodesic^2] (* Check Pythagoras *)

Here is a little bit more general code which allows to take more general curves x(t) in a surface r(u,v). We still need |x'(t)|= constant.

(* Example code to compute the geodesic curvature, O. Knill, Math 136 *) r={(3+Sin[v])Cos[u+v],(3+Sin[v])Sin[u+v],Cos[v]}; ru=D[r,u];rv=D[r,v];n=Simplify[Cross[ru,rv]];n=n/Sqrt[n.n]; x=(r/.v->s)/.u->t+s; (* here u(t)=t+s, v(t)=s for fixed s *) n=(n/.v->s)/.u->t+s; (* normal to the surface at x(t) *) velocity=D[x,t]; acceleration=D[x,{t,2}]; speed=Simplify[Sqrt[velocity.velocity]]; (* needs to be constant! *) normalcurvature = Simplify[acceleration.n/speed^2]; geodesiccurvature= Simplify[Cross[n,velocity].acceleration/speed^3]; curvature = Simplify[Sqrt[acceleration.acceleration]]/speed^2; Simplify[curvature^2-normalcurvature^2-geodesiccurvature^2]