Geodesics on a Cone

by John Oprea

This worksheet shows how to draw geodesics on the cone as well as how to calculate the arclength of a geodesic. A reference is Chapter 5 of J. Oprea, Differential Geometry and its Applications . For a=1 in , determine c and d for the geodesic connecting (1,0,1) and (0,1,1) and compare the arclength of this geodesic between these points to the arclength of the part of the parallel circle joining the points.

> with(plots):with(linalg):

Warning, the name changecoords has been redefined

Warning, the name adjoint has been redefined

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

In order to compute a geodesic, we need to know how lengths are measured on the surface of the cone. This is the metric of the surface. We can find it by computing certain tangent vectors and taking their dot products in . First, here is a dot product procedure.

> dp:=proc(X,Y)
X[1]*Y[1]+X[2]*Y[2]+X[3]*Y[3];
end:

> EFG := proc(X)

> local E,F,G,Xu,Xv;

> Xu :=[diff(X[1],u),diff(X[2],u),diff(X[3],u)];
Xv := [diff(X[1],v),diff(X[2],v),diff(X[3],v)];

> E := dp(Xu,Xu);

> F := dp(Xu,Xv);

> G := dp(Xv,Xv);

> simplify([E,F,G]);

> end:

The geodesic equations are then given by the following formulas using the metric E, F and G. In fact, the formulas below take since this is the case for the cone. The geodesic equations are a system of second order differential equations which determine a geodesic on the cone by showing how the parameters of the cone and depend on a single parameter .

> geoeq:=proc(X)

> local M,eq1,eq2;

> M:=EFG(X);

> eq1:=diff(u(t),t$2)+subs({u=u(t),v=v(t)},diff(M[1],u)/(2*M[1]))*diff(u(t),t)^2

> +subs({u=u(t),v=v(t)},diff(M[1],v)/(M[1]))*diff(u(t),t)*diff(v(t),t)

> - subs({u=u(t),v=v(t)},diff(M[3],u)/(2*M[1]))*diff(v(t),t)^2=0;

>

> eq2:=diff(v(t),t$2)-subs({u=u(t),v=v(t)},diff(M[1],v)/(2*M[3]))*diff(u(t),t)^2

> +subs({u=u(t),v=v(t)},diff(M[3],u)/(M[3]))*diff(u(t),t)*diff(v(t),t)

> + subs({u=u(t),v=v(t)},diff(M[3],v)/(2*M[3]))*diff(v(t),t)^2=0;

> eq1,eq2;

> end:

In general, the geodesic equations are impossible to solve explicitly and numerical solutions are necessary. (For the cone, an explicit solution is possible, however.) The following procedure uses Maple's numerical solution to differential equations to plot a geodesic on a surface. The inputs are: X a parametrization in u and v, ustart......vend the range of the parameters, u0,v0 the starting point of the geodesic (i.e. initial condition 1), Du0,Dv0 the starting speed of the geodesic (i.e. initial condition 2), T the range of the parameter of the geodesic itself, N a smoothness parameter describing how smooth a picture you want, gr=[d,e] where d is the number of grid lines for u and e is the number for v and orientation parameters and .

> plotgeo:=proc(X,ustart,uend,vstart,vend,u0,v0,Du0,Dv0,T,N,gr,theta,phi)

> local sys,desys,u1,v1,listp,geo,plotX;

> sys:=geoeq(X);

> desys:=dsolve({sys,u(0)=u0,v(0)=v0,D(u)(0)=Du0,D(v)(0)=Dv0},{u(t),v(t)},

> type=numeric, output=listprocedure);

> u1:=subs(desys,u(t)); v1:=subs(desys,v(t));

> geo:=spacecurve(subs(u='u1'(t),v='v1'(t),X),t=0..T, color=black,thickness=2,numpoints=N):

> plotX:=plot3d(X,u=ustart..uend,v=vstart..vend,grid=[gr[1],gr[2]],shading=XY):

> display({geo,plotX},style=wireframe,scaling=constrained,orientation=[theta,phi],shading=xy,lightmodel=light2);

> end:

Here is the cone parametrization for a=1.

> CO:=[u*cos(v),u*sin(v),u];

> EFG(CO);

> plotgeo(CO,0,1.3,0,2*Pi,1,0,-1,1,1.3,50,[5,24],45,81);

A cone geodesic condition (derived from the fact that the cone is u-Clairaut) may be found explicitly. It is
.
So, for the geodesic from (1,0,1) to (0,1,1), an angle of along the z=1 circle, and a=1, we can find D by the following. At the start and finish, u=1 and v=0 and respectively. Plugging into the geodesic formula above gives and . These can be solved to give and .

> dd:=evalf(arctan((cos(Pi/(2*sqrt(2)))-1)/sin(Pi/(2*sqrt(2)))));

> cc:=cos(dd);

So now we can plot the geodesic by replacing u in the parametrization of the cone by the formula above with the C and D just determined.

> geo:=spacecurve([cc*sec(v/sqrt(2)+dd)*cos(v),cc*sec(v/sqrt(2)+dd)*sin(v),

> cc*sec(v/sqrt(2)+dd)],v=0..Pi/2,color=black,thickness=2):

> cone:=plot3d([u*cos(v),u*sin(v),u],u=0..1.3,v=0..2*Pi,grid=[5,24],shading=xy,lightmodel=light2):

> circ:=spacecurve([cos(v),sin(v),1],v=0..Pi/2,color=red,thickness=2):

> display({geo,cone,circ},style=wireframe,scaling=constrained,orientation=[22,84]);

Here you can see geodesic in black and the parallel circle in red joining the given points. Finally, we can compute the arclength of the geodesic from (1,0,1) to (0,1,1) and compare it with the arclength of the circle parametrized by joining the two points. Of course, the circle has radius 1, so the arclength along it from (1,0,1) to (0,1,1) is simply , since the angle formed by joining the center of the circle (0,0,1) to the points is
. Here is the geodesic again.

> alpha:=[cc/cos(v/sqrt(2)+dd)*cos(v),cc/cos(v/sqrt(2)+dd)*sin(v), cc/cos(v/sqrt(2)+dd)];

We take its velocity vector and then compute the arclength integral.

> alphap:=diff(alpha,v);

> evalf(Int(sqrt((alphap[1]^2+alphap[2]^2+alphap[3]^2)),v=0..Pi/2));

So the geodesic has shorter length than the circle.