Plotting Geodesics on Surfaces

using Mathematica

Our goal is to describe in simplest terms how to use numerical integration to plot geodesics on an arbitrary surface. Geodesics are geometric invariants second in importance only to Gaussian curvature, but the differential equations for geodesics are too complicated to admit explicit solutions, except in the rarest of cases. So after some 300 years of differential geometry, geodesics were well known only for very special surfaces.

But the advent of computers has made it practical to use numerical integration to find geodesics on an arbitrary surface. We will do this with Mathematica, using its command NDSolve. Though the solutions are only approximations, they are good ones, and for the first time allow experiment with the geodesics in the most general case. (Later examples will show that precise questions about geodesics can often be solved by approximations.)

Prerequisites: Familiarity with (a) the basics of the geometry of surfaces, including geodesics, and (b) Mathematica--say the simpler commands in the Mathematica notes. In particular, for a surface in Rⁿ, 3, parametrized by x, we will need the commands ee, ff, gg that produce the metric coefficients E, F, G of the surface . Thus ee[x] evaluated on (u,v) gives the value of E at the point x[u,v].
Throughout, efg.m denotes any dot-m file (see the Mathematica notes) that contains these three commands.

The routes followed by the geodesics in a surface are strongly influenced its Gaussian curvature K. Of the three curvature commands given in the Mathematica notes, we will generally use only the simplest, gaussK, which gives the curvature of a parametrized surface x in R³. Explicitly, gaussK[x] evaluated at (u,v) gives the curvature of the surface at the point x[u,v]. We save this command in a dot-m file gaussK.m.

             Contents:
                  1. The Geodesic Differential Equations (GDE)
                  2. The Christoffel Symbols in Mathematica
                  3. Numerical Solution of the GDE
                  4. Plotting Geodesics
                  5. Examples
                  Appendix: Curvature and Geodesics

1. The Geodesic Differential Equations (GDE)

Geodesics in any geometric surface are the generalization of the straight lines in the Euclidean plane. They are the curves that do not turn. The fundamental existence and uniqueness theorem for geodesics is this:

     Given a tangent vector v at a point p of a geometric surface, there is a unique maximal geodesic that starts at p with initial velocity v. (Here "maximal" means:  of largest domain.)

We typically write (u,v) x(u,v) for the parametrization of a surface and (to avoid index ambiguity) write the square of f as f^2 rather than f². The parameter of the solutions will usually be t, with t0 its initial value.
The goal is to solve the GDE numerically for functions u(t), v(t) such that t x(u(t),v(t)) is the unique geodesic satisfying initial conditions as above. Explicitly,

• if the initial point p has x-coordinates u0, v0, then we require u(t0)=u0, v(t0)=v0, and
• if the initial velocity v has coordinates du0, dv0 (relative to the coordinate vector fields x_u, x_v), then we require u'(t0)=du0,  v'(t0)=dv0.

Thinking of (u,v) as (u₁, u₂) shows that the following geodesic differential equations (GDE) have a natural index pattern:

u'' + ¹₁₁ u^2 + 2¹₁₂ uv + ¹₂₂ v^2 = 0,

v'' + ²₁₁ u^2 + 2²₁₂ uv + ²₂₂ v^2 = 0.

The coefficients ^k_ij are called the Christoffel symbols of the parametrization. They are expressed as follows in terms of the metric components E, F, G of the parametrization, with the subscripts u and v denoting partial differentiation:

¹₁₁ = (E_uG - 2F_uF + E_vF) / (2(EG-F^2))     ²₁₁ = (2F_uE - 2E_vE - E_uF) / (2(EG-F^2))

¹₁₂ = ¹₂₁ = (E_vG - G_uF) / (2(EG-F^2))     ²₁₂ = ²₂₁ = (G_uE - E_vF) / (2(EG-F^2))

¹₂₂ = (2F_vG - G_uG - G_vF) / (2(EG-F^2))    ²₂₂ = (G_vE - 2F_vF + G_uF) / (2(EG-F^2))

Note the two equations expressing the symmetry of ^k_ij in i and j. See deriveGDE for a derivation of the geodesic differential equations and Christoffel symbols.

2. The Christoffel Symbols in Mathematica

Mathematical notation has become cryptic and sophisticated over the centuries, so the explicit unambiguous style demanded by a computer makes Mathematica's description of the Christoffel symbols a good deal longer. However, you can easily download them without typing--see below.

There are eight formulas in all:  six that appear in the GDE above and two that express the symmetry of ^k_ij in i and j.
Note: ^k_ij is given below by chr[k,i,j].

chr[1,1,1][ee_,ff_,gg_][u_,v_] := Simplify[

      (D[ee[u,v],u]*gg[u,v]-2*D[ff[u,v],u]*ff[u,v]+D[ee[u,v],v]*ff[u,v])/

      (2*(ee[u,v]*gg[u,v]-ff[u,v]^2))];

chr[1,1,2][ee_,ff_,gg_][u_,v_] := Simplify[

      (D[ee[u,v],v]*gg[u,v]-D[gg[u,v],u]*ff[u,v])/

      (2*(ee[u,v]*gg[u,v]-ff[u,v]^2))];

chr[1,2,1][ee_,ff_,gg_][u_,v_] := chr[1,1,2][ee,ff,gg][u,v];

chr[1,2,2][ee_,ff_,gg_][u_,v_] := Simplify[

      (2*D[ff[u,v],v]*gg[u,v]-D[gg[u,v],u]*gg[u,v]-D[gg[u,v],v]*ff[u,v])/

      (2*(ee[u,v]*gg[u,v]-ff[u,v]^2))];

chr[2,1,1][ee_,ff_,gg_][u_,v_] := Simplify[

      (2*D[ff[u,v],u]*ee[u,v]-D[ee[u,v],v]*ee[u,v]-D[ee[u,v],u]*ff[u,v])/

      (2*(ee[u,v]*gg[u,v]-ff[u,v]^2))];

chr[2,1,2][ee_,ff_,gg_][u_,v_] := Simplify[

      (D[gg[u,v],u]*ee[u,v]-D[ee[u,v],v]*ff[u,v])/

      (2*(ee[u,v]*gg[u,v]-ff[u,v]^2))];

chr[2,2,1][ee_,ff_,gg_][u_,v_] := chr[2,1,2][ee,ff,gg][u,v];

chr[2,2,2][ee_,ff_,gg_][u_,v_] := Simplify[

      (D[gg[u,v],v]*ee[u,v]- 2*D[ff[u,v],v]*ff[u,v]+D[gg[u,v],u]*ff[u,v])/

      (2*(ee[u,v]*gg[u,v]-ff[u,v]^2))]

The Christoffel symbols have many uses other than for geodesics, so they deserve to be saved in an individual file. The best way (though not the only way) to do this is:

• Open a new Mathematica page.
• Select--all in one batch--the formulas above.
• Copy them, and paste them into the Mathematica page. Since these formulas are separated by semicolons, the batch will be accepted as a single cell. (We gave the formulas in Plain Text rather than Bold to ease this operation.)
• Using the Cell menu, designate this cell an Initialization Cell.
• Close this page, and save it as christoffel.nb. You will then be given two choices; pick "Create Auto Save Package." This produces a file christoffel.m.

When this file is read into any later Mathematica session, by  << christoffel.m, all eight Christoffel formulas (though not displayed) will be available for use.

If for any reason, christoffel.m is not forthcoming, the standard saved version christoffel.nb will suffice. Just open it and enter the multiple command into Mathematica.

3. Numerical Solution of the GDE

Given a surface in Rⁿ, n 3, with parametrization x, we want numerical expressions for its geodesics. To get them we go through the mathematical steps that set up the GDE and arbitrary initial conditions, and for each step find a Mathematica command that will carry it out.
Then NDSolve will give the numerical solutions.

Our primary interest is in the solutions for u and v. These yield a curve t (u(t),v(t)) in the uv-plane. This curve can be moved into Rⁿ as the geodesic t x(u(t),v(t)). We call t (u(t),v(t)) the uv part of the geodesic. Obviously, the geodesic can be plotted by Mathematica only for n=3, but for n >3 the uv part is still quite informative.

      > Step 1: Get the metric components of x and name them eeh, ffh, ggh.

In the Main Case, we are given a parametrization x of a surface in Rⁿ. Now input any file efg.m that contains commands ee, ff, gg that produce the metric coefficients E, F, G of x.

Then define

       eeh[u_,v_] := ee[x][u,v];

       ffh[u_,v_] := ff[x][u,v];

       ggh[u_,v_] := gg[x][u,v]

The semicolons allow all three commands to be entered simultaneously.

The particular notation eeh, ffh, ggh is required by the command cht given below.

There is a second case where we can compute geodesics: that of an abstract geometric surface (see Section 3 of the Mathematica notes). Typically, this is a region in the plane whose metric coefficients E,F, G are given, so all we have to do is name them eeh, ffh, ggh.

      > Step 2: Substitute eeh, ffh, ggh into the general formulas for Christoffel symbols, and replace u by u[t], v by v[t].

To do this, first install the file created in Section 2 that contains the general Christoffel symbols:

     << christoffel.m

Then the required command is

      cht[i_,j_,k_][t_] := chr[i,j,k][eeh,ffh,ggh][u,v] /. {u ->u[t],v ->v[t]},

where the slash-dot operator "/." produces the substitution. (Since cht is of no independent interest, we wait and save it with a command it supports.)

      > Step 3: Set up GDE with arbitrary initial conditions, and solve numerically, using NDSolve.

A specific geodesic is determined by the GDE, the following initial conditions -- and the interval of definition of t.

• the initial value t0 of the geodesic parameter t.
• the x-coordinates u0,v0 of the initial point p, so p=(t0)=x(u0,v0),
• the x-coordinates du0,dv0 of the initial velocity v, so '(t0)= du0 x_u(u0,v0) + dv0 x_v(u0,v0),
• the end points tmin and tmax of the interval tmin < t < tmax on which we seek solutions (this interval must contain t0).

Now we have some good luck: Because the functions u(t) and v(t) do not appear explicitly in the GDE, its two second-order DEs can be recast more simply as four first-order DEs. It suffices to take the first derivatives u' and v' as new variables: p=u', q=v'. Then obviously, p'=u", q'=v", and the GDE take this form:

                                 u'=p

                                 v'=q

                                 p' + ¹₁₁ u^2 + 2¹₁₂ uv + ¹₂₂ v^2 = 0

                                 q' + ²₁₁ u^2 + 2²₁₂ uv + ²₂₂ v^2 =0

We have seen how to write these equations in terms of Mathematica and thus the following command gdesolv will produce a numerical solution of the GDE satisfying the designated initial conditions.

    gdesolv[u0_,v0_,du0_,dv0_,t0_,tmin_, tmax_] :=

         NDSolve[ {u'[t]==p[t], v'[t]==q[t],

     p'[t]+cht[1,1,1][t]*p[t]^2+2cht[1,1,2][t]*p[t]*q[t]+

          cht[1,2,2][t]*q[t]^2==0,

     q'[t]+cht[2,1,1][t]*p[t]^2+2cht[2,1,2][t]*p[t]*q[t]+

          cht[2,2,2][t]*q[t]^2==0,

    u[t0]==u0, v[t0]==v0, p[t0]== du0, q[t0]==dv0},

         {u,v,p,q}, {t,tmin,tmax}].

Note here that the differential equations and the initial conditions are all lumped together within a single pair of braces {...}. As usual, double equal signs are used for equalities that are sought rather than known.

The solutions produced by this command are not direct functions u,v,p,q of t, but rather Interpolating Functions--data that produces requested values of u,v,p,q by interpolation.

To avoid confusion it's a good idea to assign names to the outputs of gdesolv, say, soln1, soln2,... or others more descriptive.

Now we save this command gdesolv and the command cht from the preceding section, together in a one file: gdesolv.m, which is constructed just as for christoffel.m, but with two initialization cells.

4. Plotting the Geodesic

A numerical solution from Step 3 above, say soln4, is plotted in the uv-plane by

ParametricPlot[Evaluate[{u[t],v[t]}/.soln4], {t,tmin,tmax}]

"Evaluate[...]" speeds things up and prevents Mathematica from complaining. To get a better picture, it will sometimes be necessary to rerun gdesolv to cover a larger interval for t.

Various options are useful in plotting.

• For plots in dimension 2, inserting AspectRatio ->Automatic at the end of the command above forces Mathematica to use the same scale on both axes.
• ImageSize ->pqr for numbers such as 216 or 360, where 72=1 inch. (The default size is 288.)
• For more choices, such as PlotRange, type: Options[ParametricPlot].

In the Main Case, for all n 3, we get an InterpolatingFunction solution for the geodesic t x(u(t),v(t)) in Rⁿ. When n=3, Mathematica can plot it, using:

     ParametricPlot3D[ Evaluate[x[u[t],v[t]] /. soln4], {t,tmin,tmax}]

In the Abstract Case, the curve t (u(t),v(t)) in R² is already the required geodesic.

5. Two Examples

We test the commands derived above in the Main and Abstract cases..

The Main Case
Let's find a geodesic on the unit sphere S in R³. The geodesics on the sphere follow great circles, so geometrically speaking, there is only one.

We parametrize S by x: DR³, where x(u,v) = (cos(v) cos(u), cos(v) sin(u), sin(v)) and the region D is given by - u ,   -/2 v /2.

Now, in a new Mathematica page, type the Mathematica version of x and install a file that efg.m will compute the metric components of x.

     x[u_,v_] := {Cos[v]*Cos[u], Cos[v]*Sin[u], Sin[v]}

     << efg.m

Folowing the scheme in the previous section, we first define

     eeh[u_,v_] := ee[x][u,v];   ffh[u_,v_] := ff[x][u,v];   ggh[u_,v_] := gg[x][u,v]

and then install the two dot-m files constructed there.

     << christoffel.m

     << gdesolv.m

We can now solve the GDE for a geodesic of S:

     soln1= gdesolv[0,0,1/2,Sqrt[3]/2,0,-4,4]

where the numbers 0,0 locate the initial point of the geodesic at x(0,0) and the next two determine its initial velocity. The final three require that (0)=p and that is defined on the interval -4 < t < 4. (This interval should cover what will actually be needed.)

[This is the first command that produces a visible response, a few lines about InterpolatingFunctions. Recall that in the general definition of gdesolv, it is evaluated on u0_,v0_,du0_,dv0_,t0_,tmin_,tmax_. The coordinate vectors x_u, x_u are orthonormal at (u,v)=(0,0), so we have asked for a unit-speed geodesic starting at x(0,0).]

ParametricPlot[Evaluate[{u[t],v[t]}/.soln1], {t,-Pi,Pi}, AspectRatio -> Automatic, ImageSize ->288]

      [This should produce what looks like a sin curve. The values -Pi,Pi are clear since a great circle on S has length 2.]

ParametricPlot3D[Evaluate[x[u[t],v[t]]/.soln1], {t,-Pi,Pi}, ImageSize ->252]

      [If all has gone well, this will give :]

The Abstract Case
Consider the Poincaré disk model of the hyperbolic plane (Section 3 of the Mathematica notes). This an abstract geometric surface defined on the disk D: r< 2 in R² with inner product characterized by metric components E, F, G such that

       eeh[u_,v_] := 1/(1-(u^2+v^2)/4)^2

       ffh[u_,v_] := 0

       ggh[u_,v_] := 1/(1-(u^2+v^2)/4)^2

Here the metric components are already named eeh,ffh,ggh.

For next two steps, load christoffel.m and gdesolv.m (which also contains the command cht).

For the geodesic starting at, say, (u,v)=(1,0) in the direction (0,1), the general format

      gdesolv[u0_,v0_,du0_,dv0_,t0_,tmin_,tmax_]

is filled in by

     soln1= gdesolv[1,0,0,1,0,-20,20]

Another geodesic is computed by

     soln2= gdesolv[1,0,-1,0.5,0,-20,20]

Now plot these by

    g1=ParametricPlot[Evaluate[{u[t],v[t]}/.soln1],{t,-20,20}],

then g2 similarly. For reference, we also want the rim of the disk.

    rim= ParametricPlot[{2*Cos[t],2*Sin[t]},{t,0,2Pi}]

These figures will be misshapen, but that can be corrected in the combined picture:

    Show[{g1,g2,rim}, AspectRatio ->Automatic, ImageSize ->180]

In fact, it is known that the geodesics of this surface follow circular arcs orthogonal to the rim, with straight lines through the origin as a limiting case, so the figure shows four geodesics.

Appendix:  Curvature and Geodesics

The two most important geometric invariants, Gaussian curvature and geodesics, are related in various ways. Here is one of the most direct; it will supply remarkable information about where geodesics go and why.

In an oriented surface (such as R²), let D be a simply connected region whose boundary D consists of a finite number of geodesic segments. Consider one trip around this boundary in the positively oriented direction, that is, "keeping D on the left" in the R² case.

At each vertex--where two segments meet--the turning angle is the smallest-in-absolute-value angle , measured in radians, through which the velocity vector of the incoming segment must turn to reach that of the outgoing segment.
The figure below shows the usual case where the turn is in the positive direction, so >0. (For the negative direction, <0.)

Then the so-called "Gauss-Bonnet theorem with angles" ([EDG],[G]) gives:

      Corollary.    TK(D) + TA(D) = 2,

where TK(D) is the total curvature of D, that is, the integral of Gaussian curvature K over D, and
TA(D) is the total turning angle, that is, the algebraic sum of all the turning angles on the boundary of D.

Examples. (1) If D is a hemisphere in a sphere S of radius r, its boundary is a (geodesic) great circle. Thus there are no turning angles, and the formula predicts TK(D)=2. In fact, since S has constant curvature K=1/r²,

TK(D)=(1/r²) area(D)= (1/r²)(4r²/2)= 2.

(2) In the figure at left, D is the interior of a geodesic loop in an abstract geometric surface. At its single vertex, the turning angle is the angle from final velocity w of the boundary geodesic to its initial velocity v.

Without knowing the metric components of this surface, we cannot compute this angle. But we can be sure that it is less than . The Corollary above then shows that TK(D) > . Thus such geodesic loops can exist only where there is considerable positive curvature.

(3) If a surface has K0, then two geodesics starting from the same point cannot meet again to enclose a simply connected region D. In fact, such a region must have TK(D)>0, again since both turning angles on D are less than .

HOME