Geodesics on Two Hills

Geodesics on a Surface with Two Hills

Our goal is to plot the geodesics of a naturally chosen one-parameter family of geodesics on a surface with two unequal hills. We follow the method described in The Setup.

           Contents
               1. The Hills
               2. Gaussian Curvature
               3.The Plotting Machine
               4. The Critical Geodesic G1
               5. Total Curvatures
               6.The Critical Geodesic G2
               7. A Sample Geodesic in 3-space
               8. The Critical Geodesic G3 (compare G2)
               9. The Critical Geodesic G4 (compare G1)
              10.The Critical Geodesics G5 and G6
              11. Summary

1. The Hills

The surface is parametrized by these formulas:

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

x[u_,v_]:= {u,v,h[u,v]}

This is a Monge parametrization, completely described by the function h that gives the height of each point above the uv-plane. The surface is shown by:

hills = ParametricPlot3D[Evaluate[x[u,v]],{u,-4,4.6},{v,-2.5,2.5},

PlotRange->All, PlotPoints->40, Boxed->False]

The taller hill, called hill1, reaches a height of about 4, near (2,0). The smaller hill, hill2, has height about 2, near (-2,0). Here, as is often the case, there is no need for more precision.

There is a natural geodesic in this surface, namely, the curve that lies over the u-axis By symmetry, it is a geodesic; we call it the profile geodesic. It is plotted, though not yet shown, by the command:

profilegeod = ParametricPlot3D[Append[{t,0,h[t,0]},

AbsoluteThickness[2]]//Evaluate,{t,-4,4.6},Boxed->False,

Axes->False, DisplayFunction->Identity]

Then it is shown on the surface by:

Show[{hills,profilegeod}, PlotRange->All, ViewPoint->{1.3,-1.4,2},

Axes->False, DisplayFunction->$DisplayFunction]

We can now state our goal more precisely: it is to organize and understand the one-parameter family of all geodesics orthogonal to the profile geodesic. The parameter that distinguishes among these geodesics is the value of u at which they cross the profile curve.

Probably not many people could guess, ahead of time, what these geodesics will look like, simply because, until recently, there was no reasonable way to determine where a particular geodesic would go.

2. Gaussian Curvature

To compute Gaussian curvature, we install the command gaussK from the Mathematica notes, but any equivalent will do as well.

<< gaussK.m

The following operations produce a function kkf whose value at (u,v) is the Gaussian curvature K of the surface at the point x(u,v).

kkuuvv = gaussK[x][uu,vv];

kkf[u_,v_]:= kkuuvv /.{uu->u,vv->v}

The following race will show the advantage of kkf:

gaussK[x][1.1,-0.5]//Timing

kkf[1.1,-0.5]//Timing

A fundamental fact about this surface is that as distance from the hills increases, it becomes ever flatter. The surface is said to be "flat at infinity." One way to see this is to plot K on a few large circles. For example, at radius 100,

Plot[kkf[100*Cos[t],100*Sin[t]], {t,0,2Pi}]

gives values of order 10^-14. (Try radius 10⁶.)

Here is a key picture that shows the relationship of height h to curvature K along the profile geodesic.

Plot[{h[u,0],kkf[u,0]}//Evaluate,{u,-4,4.6},

PlotStyle-> {RGBColor[0,0,1],RGBColor[1,0,0]}]

The (red) curvature graph is cut off at the top since K is so large on the hilltops. In fact, N[kkf[2,0]] shows that K is above 63 near the top of hill1, and similarly K is over 15 at the top of hill2.

The contour plot below, roughly aligned with the graphs, indicates the distribution of curvature away from the profile geodesic.

The figure-8 region has the lowest values of K<0. The white disks-- the tops of the hills--are the only regions with K>0, except for some microscopically small values along the v-axis.
But on the hilltops, the K>0 curvature is very large, while nowhere is the K<0 curvature substantial. (See the red graph above.)

As will soon appear, the isolated, highly curved hilltops have a powerful effect on the geodesics of the surface. The simplest way see how the the hilltops are related to the paths followed by geodesics is to draw the two circles that roughly outline them.

hilltops = ParametricPlot[{{2+.576*Cos[t],.576*Sin[t]},

{-2+.57*Cos[t],.57*Sin[t]}}//Evaluate,{t,0,2Pi},

PlotStyle->{{RGBColor[0,0,1],AbsoluteThickness[2]},

{RGBColor[0,0,1],AbsoluteThickness[2]}},AspectRatio->Automatic]

The two hilltops lie over disks of about the same size in the uv-plane, but recall that the top of the taller hill1--on the right--has curvature much greater than that of hill2.

3. The Plotting Machine

We apply the approach in The Setup to get geodesics orthogonal to the profile geodesic.

The first step is to input a file such as efg.m from the Mathematica notes that will compute the metric coefficients of the parametrization x.

<< efg.m

Name these as follows (necessary to fit later commands):

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

Next we input the commands needed for geodesics.

<< christoffel.m

<< gdesolv.m

Recall that the command gdesolv feeds on initial conditions and parameter data as follows:

gdesolv[u0,v0,du0,dv0,t0,tmin,tmax]

To produce geodesics orthogonal to the profile geodesic, we take

(u0,v0) = (u0,0), (du0,dv0) = (0,1), t0 = 0, tmin = -b, tmax = b.

Here is a machine that will numerically integrate these geodesics:

mac[u0_,b_]:= gdesolv[u0,0,0,1,0,-b,b]

Entering ggh[u,0] yields 1, which implies that our geodesics initially have unit speed. But geodesics have constant speed, so they have unit speed everywhere, that is, they have arclength parametrization.

4. The Critical Geodesic G1

To move through the one-parameter family of geodesics, we elect to have the distinguishing parameter u0 of our geodesics start with large positive values and decrease. For u0 large, the geodesic through (u0,0) is almost a vertical line since K is already so close to 0 that the geometry is almost Euclidean.

As we move u0 downward toward hill1, the large positive curvature there takes effect on the routes followed by the geodesic. To illustrate this we plot the geodesics with four successive values for u0: 3.5, 3.25, 3.0, 2.75.

The computations for these geodesics can be combined into a single command by separating the individuals by semicolons, as follows:

soln3p5 = mac[3.5,8]; soln3p25 = mac[3.25,8];

soln3p0 = mac[3.0,8]; soln2p75 = mac[2.75,8];

Then they are plotted all at once (but not yet shown) by:

fourgeods = ParametricPlot[{{u[t],v[t]}/.soln3p5,{u[t],v[t]}/.soln3p25,

{u[t],v[t]}/.soln3p0,{u[t],v[t]}/.soln2p75},

{t,-8,8},DisplayFunction->Identity]

The use of "Evaluate" in the preceding command may speed the computation. Now we can show the geodesics--along with the hilltop circles--by the command:

Show[{fourgeods,hilltops}, AspectRatio->Automatic, PlotRange->All]

Here we have cropped the picture at v=±4.25. No information is lost since as soon as a geodesic gets away from the hills, it becomes virtually straight (again, because curvature has become negligibly small).

Reading these geodesics from right to left suggests that something new is going to happen as u0 continues descending past 2.75.

The geodesic with u0=2.71 is not radically different, but the geodesic with u0=2.70 is a different animal than those that came before: they were all one-to-one, but this one cuts across itself near u=-25. We show these last two geodesics (below) without using the option "AspectRatio->Automatic". This frees Mathematica to use different scales on the axes, and in many cases this gives more informative pictures.

soln2p71 = mac[2.71,35];soln2p70=mac[2.70,35];

ParametricPlot[ {{u[t],v[t]}/.soln2p71,{u[t],v[t]}/.soln2p70},

{t,-35,35}, PlotRange->All, ImageSize->216]

Choosing values of u0 back and forth in the interval 2.70< u0 <2.71, we find analogous pairs of geodesics on ever shorter sub-intervals. These suggest correctly that there is a smallest number u0[1] whose geodesic is still one-to-one, while all those with smaller u0 have crossings.

A good approximation of u0[1] turns out to be 2.7073, so we plot this geodesic:

G1 = ParametricPlot[Evaluate[{u[t],v[t]}/.soln2p7073],

{t,-25,25}, PlotRange->All]

The branches here show as parallel from -5 out to -30, and since the Gaussian curvature at (-30,0) is already of order 10^-10, they are likely to remain parallel on out to -oo. But even the small drop of 0.0001 to u0=2.7072 yields a geodesic that has a self-intersection. Now 2.7073 is doubtless not the exact cut value, but it suffices to signal a change in character of the geodesics.

The geodesic G1, with u0[1] = ^~2.7073, is the first critical geodesic. G1 is an infinite geodesic loop, with vertex at u=-oo. It encloses a simply connected region D that includes a long tail of width ^~2.5 out to -oo.

5. Total Curvatures

We consider some consequences of the Gauss-Bonnet corollary in the Appendix of The Setup. The critical geodesic G1 can be regarded as the limit of geodesics with u0 approaching u0[1] from below. (See g2p70 in the comparison figure above.) Thus G1 can legitimately be considered to have turning angle

at -oo. So the Gauss-Bonnet relation for the region D1 is TK +

= 2

, giving TK=

for the total curvature inside the loop.

Is this reasonable? Well, let Mathematica estimate the total curvature of the hilltops. Using x[2+r*Cos[q],r*Sin[q]] on hilltop 1, numerical integration of K dA gives ^~4 for its total curvature. Similarly, we get ^~2.4 for hilltop 2. Thus in the preceding figure, a large area of (small) K<0 curvature is needed in D1 to reduce the total positive curvature ^~6.4 to the net of .

Incidentally, the entire surface has total curvature 0. This can be shown by applying the Gauss-Bonnet corollary to increasingly large near-squares, whose four turning angles will each approach the Euclidean value /2. Thus the vast region outside G1 has mostly negative curvature, but so small that its total is only -.

6. The Critical Geodesic G2

We continue the descent of u0 from 2.7 down to 2.6. There

soln2p6=mac[2.6,8]; g2p6=ParametricPlot[Evaluate[

{u[t],v[t]}/.soln2p6], {t,-8,8},DisplayFunction->Identity]

leads to

Show[{g2p6,hilltops},AspectRatio->Automatic,ImageSize->252,

DisplayFunction ->$DisplayFunction]

The crossing point, near -25 for g2p7 has moved in to near -2.4, on the slope of hill2. A fair guess is a critical geodesic may occur when u0 is slightly lower, making the geodesic tangent to the top circle of hill2, say at u0=2.57. But this value is already too low.

Successive over and under estimates, as in the search for G1, finally yield

soln2p57464 = mac[2.57464,30]; G2 = ParametricPlot[

Evaluate[{u[t],v[t]}/.soln2p 57464],{t,-30,30},PlotRange->All]

We accept this geodesic (above left) as critical since, as before, the two branches show as parallel out to negligibly small values of curvature. Above right is an enlargment of the region near the crossing of the u-axis.

G2 will look different when, as below), equal scales are imposed on both axes and the hilltops are included.

The geodesic G2, with u0[2]=^~2.57464, is the second critical geodesic. It has a crossing near the top of hill2, and then a long thin strip that runs out to a limit vertex at u=-oo. G2 encloses two simply connected regions: Dh, which includes hill1, and Ds, the strip. Ds is so much thinner than the analogous strip for G1 that it doesn't show up well on the same-scale figure below.

Show[{G2,hilltops}, AspectRatio->Automatic,

PlotRange->All, ImageSize->360]

To find the total curvature of these regions we need their turning angles at their shared vertex (the crossing point). With the regions oriented in the usual way, this turning angle is the same for both. A closeup view of the angle (with "AspectRatio->Automatic" in force) suggests an angle moderately less than radians, probably not too far from 3. We compute it to be ^~2.973.

By the Gauss-Bonnet corollary:

TK(Dh)+2.973= 2, so TK(Dh)=^~3.310,

Plausible, since hilltop1 has TK approximately 4.

TK(Ds)+2.973+ =2, so TK(Ds)= -2.973 = ^~0.1685.

This small positive number is not unreasonable since the bulge of Ds is on hilltop2, where K>0 is large, while on the long tail, K<0 is extremely small.

Note that TK(Dh) = +TK(Ds).

7. A Sample Geodesic in 3-Space

Any reduction of u0 below u0[2]=2.57464 will pull the vertex at u=- in to a finite value. By u0=2.57462 the vertex is all the way in to about -12.5.

As a sample of this kind of geodesic, take u0=2.55

soln2p55=mac[2.55,20]; g2p55=ParametricPlot[Evaluate[

{u[t],v[t]}/.soln2p55], {t,-8,8},DisplayFunction->Identity]

The geodesic is plotted, with the hilltops for reference, by:

Show[{g2p55,hilltops},AspectRatio->Automatic,

DisplayFunction->$DisplayFunction]

Note that the second crossing of the u-axis has closed in to about 2.5.

We want to show this geodesic up on the surface in 3-space. Mathematica does not handle this assignment particularly well if either the surface or the geodesic is at all complicated. The trouble comes from the way a surface is plotted: as a collection of flat opaque panels. A geodesic can be plotted more accurately, but any of its pixels behind a panel disappear from view.

As an alternative to the standard plot of our surface, we will make a grid model consisting of a few u-parameter and v-parameter curves. The intervening space between these curves is not filled in--it's just a grid.

The sample geodesic is lifted into R³ and its line thickened by the command:

sam=ParametricPlot3D[Evaluate[Append[x[u[t],v[t]],

AbsoluteThickness[3]] /.soln2p55], {t,-10,10}]

Now we drape the geodesic over the grid and try for a good viewpoint.

Show[{grid,sam},AspectRatio->Automatic, Boxed->False,Axes->False,

PlotRange->All, ViewPoint->{-.4,-.8,.8}]

Replace "grid" by "hills" to get the second figure below.

Note: Figures within Mathematica look and print better than those exported to the internet, so we have had to slightly edit these and other figures.

8. The Critical Geodesic G3

Descending from u0=2.55, try:

soln2p2 = mac[2.2,8]; ParametricPlot[Evaluate

[{u[t],v[t]}/.soln2p2], {t,-8,8},AspectRatio->Automatic]

It might seem as if the first vertex (near 0) is going to be pulled apart by lifting the v>0 curves and lowering the v<0 curves, but that is impossible since what we have here is a crossing. As u0 continues to decline, the jaws opening toward u=-oo gradually close and a critical geodesic is reached that resembles G2 (see below).

The geodesic G3, with u0[3]=^~2.19128, is the third critical geodesic. Like G2, it has a crossing near the top of hill2, and then a long thin strip out to a limit vertex at u=-oo. And like G2, it encloses two simply connected regions: Dh, which includes the top of hill1, and the strip Ds.

Both G3 and G2 surround the same total curvature, but G3 encloses less of hill1than G2, hence it must have a smaller K<0 area. Both have strips out to u=-oo that are too thin to show up well in same-scale figures.

9. The Critical Geodesic G4

With only a very small decrease from u0[3], the edges of the infinite strip open back up and the crossing moves to lower values of u, so that g2p15 looks like this:

As u0 contines to descend, the crossing goes out to infinity, yielding another critical geodesic:

The geodesic G4, with u0[2]=^~2.19128, is the fourth critical geodesic. Like G1, it is a geodesic loop with a limit vertex at u=-oo. It encloses hill2 as does G1, but unlike G1 (see below), it encloses only part of hill1. As with all such geodesic loops, the total curvature is .

We show both G1 and G4, along with the hills by the command:

Show[{G1,G4,hilltops},AspectRatio->Automatic,

PlotRange->{{-5,3},{-1.4,1.4}}]

Note that since G1 and G4 each surround a region of total curvature , the region between them has total curvature 0.

10. The Critical Geodesics G5 and G6

As soon as u0 drops below u0[4], parallel lines of G4 open to give the same kind of curves as for u0>u0[1],that is, curves shaped like a parabola opening to the left. They straighten out as u0 approaches the top of hill1 and then become parabola-shaped opening to the right.

Up until now, all the geodesics have moved in the -u direction as they left their initial points. But here as the initial point crosses the top of hill1 at u0=^~1.99, the orientation reverses. Once past the top of hill1, its "gravitational pull" tugs geodesics back in the +u direction.

As u0 moves on down the slope of hill1,the roughly parabolic shapes squeeze down as in earlier cases to a critical geodesic

soln1p81523=mac[1.81523,20]; G5=ParametricPlot[Evaluate[

{u[t],v[t]} /.soln1p181523],{t,-20,20}, PlotRange->All]

This passes the parallelism test, so is a critical geodesic.

The geodesic G5, with u0[5] = ^~1.81523, is the fifth critical geodesic. It is again an infinite geodesic loop, but its vertex is at +oo, thus reversing the infinite loops G1 and G4. Like all infinite loops, it encloses a region of total curvature .

Show[{G5,hilltops},PlotRange->All, AspectRatio->Automatic]

Values of u0 immediately below u0[5] have self-intersections on the positive u-axis.

As u0 decreases, the crossing point moves closer to u=0, with a closest approach somewhere near 3 when u0 is near 1.7, as in the picture to the right.

As u0 continues to descend, the crossing point moves back out again until it reaches +oo producing another critical point.

The geodesic G6, u0[6]=^~1.2786, is the sixth critical geodesic. Like G5, it is an infinite geodesic loop with vertex at u=+oo and encloses a region of total curvature .

G6 has descended the slope of hill1 until it encloses the entire top of hill1. Since G5 and G6 enclose regions of the same total curvature, G6 will require more area of negative curvature than G5. This is supplied by the width of its strip out to +oo, which at ^~1.4 is four times that of G5.

We have found six critical geodesics that start on hill1, so now the question is: How many will hill2 yield?

The answer is none. Its smaller total curvature (2.4 vs 4.0 for hill1) can deflect geodesics with u0<u0[6], but evidently is not strong enough to produce self-intersections.

The following rather crowded picture shows a sequence of geodesics with u0 descending from near u0[6] down to below -3.

As u0 progresses toward hill2, the only noteworthy change is that geodesics initially opening to the left switch to the right, around the minimum point u=-0.3 of h[u,0].

Then as u0 passes the hilltop near -2, the geodesics quickly reverse direction again (red curve). Further down at u0=2.56, the geodesic is at the edge of the K>0 hilltop, where presumably the bending power of hill2 is strongest. But no self-intersections result. Then as u0 moves away from hill2, the geodesics flatten out, since curvature is rapidly approaching 0.

11. Summary

A natural way to classify geodesics is by the number # of their self- intersections. In our case--the family of all geodesics starting orthogonal to the profile geodesic--these numbers are 0, 1, 2. The initial u-coordinates u0[1],...,u0[6] of six critical geodesics mark off the intervals where # is 0, 1, or 2.

This information is summarized in the graph to the right, where # is plotted against u0. The height function uh[u,0], shown in blue, outlines hill1 (not to scale).

The graph shows that all the action starts on hill1. As u0 varies, the number #[u0] only has jumps of ± 1: there is no way for one-to-one geodesic to suddenly acquire two intersections, or the reverse.

The critical geodesics occur in pairs:
   G1 and G4 are geodesic loops with a vertex at u=-oo,
   G5 and G6 are similar but reversed, so their vertices are at u=+oo,
   G2 and G3 enclose two regions and have one finite vertex and one at u=-oo.

We could further distinguish one-to-one geodesics that initially move in the +u dir from those moving in the -u direction, for these transitions occur only at the three critical points of the height function on the profile geodesic.

HOME