# How do you find a geodesic on an ellipsoid?

This question has been asked twice but not been answered on this site, and I really need to understand this concept.

All I want is, based on 2 arbitrary points A and B on a given ellipsoid, to find the geodesic containing those 2 points.

• Since some ellipsoids are surfaces of revoution, en.wikipedia.org/wiki/Clairaut%27s_relation helps really much in such cases. Commented Oct 22, 2017 at 2:08
• @Jack D'Aurizio: two of the axes should be equal for that; but it's a good start, past the sphere case. Commented Oct 22, 2017 at 2:15
• I wrote a full-on honors research paper on this freshman year. I'll go through it now and write as in-depth of an answer as I can, please don't delete this question. Commented Oct 22, 2017 at 2:16
• Andrew Tawfeek, thank you. I won't. Commented Oct 22, 2017 at 2:20
• It looks like the first solution to the triaxial case came from Jacobi: en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid Commented Oct 22, 2017 at 2:32

So a geodesic is essentially the minimization of a functional, that is, a function of functions. They are the main study of Calculus of Variations. Now, if we have a functional of the form $$J[y]=\int_a^b F(x,y,y') \ dx$$ defined on the set of functions $y(x)$ which have continuous first derivatives in $[a,b]$, and satisfy the boundary conditions $y(a)=A$, $y(b)=B$, then a neccesary condition for $J[y]$ to have an extremum for a function $y(x)$ is that $y(x)$ satisfy the Euler-Lagrange equation $$F_y-\frac{d}{dx}F_{y'}=0.$$

In essence, we need to solve that above differential equation.

Without getting anymore nitty gritty on things, the functional we care about here is the arc length integral.

Firstly, suppose we have a surface $\mathcal{S}$ defined by the vector equation $$\vec{r}(u,v)=x(u,v) \hat{i} + y(u,v) \hat{j} + z(u,v) \hat{k}.$$ The geodesic curve lying on surface $\mathcal{S}$ can be specified by the equations $$u=u(t), \qquad v=v(t),$$ and can be found by minimizing the arc length integeral $$L=\int_{t_0}^{t_1} \sqrt{\left( \frac{dx}{dt} \right)^2+\left( \frac{dy}{dt} \right)^2+\left( \frac{dz}{dt} \right)^2} \ dt$$ But since $x,y,z$ are each functions of more than one variable, the chain rule is going to turn this integral into a monstrosity. Thankfully, it's a monstrosity with a pretty neat pattern so we can make things look a bit prettier.

Letting $u'=\frac{du}{dt}$ and $v'=\frac{dv}{dt}$, the integeral above can just be rewritten as

$$J[u,v]=\int_{t_0}^{t_1} \sqrt{Eu'^2+2Fu'v'+Gv'^2} \ dt,$$

where $E$, $F$, and $G$ are the coefficients of the first fundamental (quadradic) form of the surface, i.e., $$\begin{cases} \displaystyle{E=\vec{r}_u\cdot \vec{r}_u = \left( \frac{\partial x}{\partial u} \right) ^2 +\left( \frac{\partial y}{\partial u} \right) ^2 + \left( \frac{\partial z}{\partial u} \right) ^2} \\ \\ \displaystyle{F=\vec{r}_u \cdot \vec{r}_v = \frac{\partial x}{\partial u} \frac{\partial x}{\partial v} + \frac{\partial y}{\partial u} \frac{\partial y}{\partial v} + \frac{\partial z}{\partial u} \frac{\partial z}{\partial v}}\\ \\ \displaystyle{G=\vec{r}_v\cdot \vec{r}_v=\left( \frac{\partial x}{\partial v} \right) ^2 +\left( \frac{\partial y}{\partial v} \right) ^2 + \left( \frac{\partial z}{\partial v} \right) ^2}. \end{cases}$$

The Euler-Lagrange equation in this case corresponds to the two different equations $$\displaystyle{ F_u- \frac{d}{dt} F_{u'} = 0,} \ \ \ \ \ \ \ \displaystyle{F_v- \frac{d}{dt} F_{v'} = 0,}$$ hence, we obtain $$\label{geodesiceq1} \frac{E_u u'^2+2F_u u' v' + G_u v'^2}{\sqrt{Eu'^2+2Fu'v'+Gv'^2}} - \frac{d}{dt} \frac{2(Eu'+Fv')}{\sqrt{Eu'^2+2Fu'v'+Gv'^2}} = 0,$$ $$\label{geodesiceq2} \frac{E_v u'^2+2F_v u' v' + G_v v'^2}{\sqrt{Eu'^2+2Fu'v'+Gv'^2}} - \frac{d}{dt} \frac{2(Fu'+Gv')}{\sqrt{Eu'^2+2Fu'v'+Gv'^2}} = 0.$$

Which are the two differential equations whose solutions provide the geodesic on surface $\mathcal{S}$.

Now let me note one important case. The case where $E$ and $G$ are explicit functions of $v$ only and $F=0$, we have $$\frac{E_v+v'^2G_v}{2\sqrt{E+Gv'^2}}-\frac{d}{du} \left( \frac{Gv'}{\sqrt{E+Gv'^2}} \right) =0,$$ so $$E_v + v'^2G_v - 2 G \sqrt{E+Gv'^2} \bigg[ \frac{v''}{\sqrt{E+Gv'^2}} + \left( \frac{1}{2} \right) \frac{v'(2Gv'v'')}{(E+Gv'^2)^{\frac{3}{2}}} \bigg] =0$$ $$E_v + v'^2 G_v - 2Gv'' + \frac{2G^2v'^2v''}{E+Gv'^2} = 0$$ $$\frac{Gv'^2}{\sqrt{E+Gv'^2}} - \sqrt{E+Gv'^2} = c_1$$ which can be made even more helpful by noting $\displaystyle{v'= \frac{dv}{du}}$, giving us $$Gv'^2-(E+Gv'^2)=c_1 \sqrt{E+Gv'^2}$$ $$\left( - \frac{E}{c_1} \right) ^2 = E + Gv'^2$$ $$\frac{E^2-{c_1}^2E}{G{c_1}^2} = v'^2,$$ finally providing $$u={c_1} \int_{v_1=v(u_1)}^{v_2=v(u_2)} \sqrt{ \frac{G}{E^2-{c_1}^2E }} \ dv.$$

Now lets talk about the ellipsoid. An ellipsoid with semi-axis lengths $a,b,c$ is given by the vector equation $$\vec{r}(u,v)= \langle a \cos(u) \sin(v), b \sin(u) \sin(v), c \cos(v) \rangle$$ with $u \in [0,2\pi)$ and $v\in [0, \pi]$.

The first fundamental form of the ellipsoid which I sadly found by hand before I realized I could've went to Google for is

$$\begin{cases} E= b^2 \cos^2(u) \sin^2(v)+a^2 \sin^2(u) \sin^2(v) \\ F= \left( b^2-a^2 \right) \sin(v) \cos(v) \sin(u) \cos(u)\\ G= \cos^2(v) \left[ (a^2-b^2) \cos^2(u)+b^2 \right] +c^2 \sin^2(v) \end{cases}$$

Now listen, I'm actually procrasinating the hell out of a transfer essay by being here and looking at the fundamental form above is pretty worrying, because it's far from elegant and being easy to deal with and would entail solving that frightening differential equation I mentioned earlier.

So I'm going to let $a=b$. Playing around with $c$ is ellipsoid-y enough for me, and if you want to go for the full $a,b,c$ then be my guest, haha! So doing what I said, things are a lot prettier:

$$\begin{cases} E= a^2 \sin^2(v) \\ F= 0 \\ G=b^2 \cos^2(v) + c^2 \sin^2(v) \end{cases}$$

and this is where we can apply that last remark I made earlier regarding the first fundamental form looking like this. We thus have the geodesic on this surface given by

$$u={c_1} \int \sqrt{ \frac{G}{E^2-{c_1}^2E }} \ dv = c_1 \int \sqrt{\frac{b^2 \cos^2(v) + c^2 \sin^2(v)}{a^4 \sin^4(v) - c_1^2 a^2 \sin^2(v)}} \ dv.$$

Now in the case of a sphere, solving the integeral and rearranging gives a plane -- and that planes interesection with the sphere represents geodesics on a sphere. I myself couldn't get anywhere with integeral, but I'm sure someone else can so I'll just leave things off here.

... Now back to my essay.

On a very final note, here is my research paper from my first semester. I prove every step up to proving the very Euler-Lagrange equations themselves and then show some examples of it in use, like finding the geodesic on a sphere and geodesic of a function revolved around an axis.

Like Jack mentioned in the comments to your post, if you consider just the ellipsoids that can be obtained through revolving a function $f(x)$, the integeral is much easier. I prove it in the above paper, but here it is again:

$$v=c_1 \int_{u_1}^{u_2} \frac{\sqrt{1+(f'(u))^2}}{f(u) \sqrt{(f(u))^2-{c_1}^2}} \ du$$

Note, $x=u$ in the above integeral.

• From your calculation, your $G$ should be (after setting $a=b$) $G=b^2\cos^2(v)+c^2\sin^2(v)$?
– user99914
Commented Oct 22, 2017 at 4:38
• @JohnMa Ahh yes, you're correct! Sorry about that, let me edit it now. Commented Oct 22, 2017 at 4:48
• Thank you so much, I am an undergrad in computer science, and this is hard to swallow right now so i am just going to stare at it for a while doing my best to understand it. Thank you very much. As a side question, please tell me you are at least a Grad student so that my ego doesn't plummet to the ground :p Commented Oct 22, 2017 at 4:55
• @Mak, I will say that explaining elliptic integrals to a five-year-old looks like an unreasonable request. ;) Nevertheless: these are nonelementary functions that will inevitably show up when you're integrating something that came out of an ellipse or an ellipsoid. Most modern computing environments should have a stable implementation of them. Commented Oct 22, 2017 at 6:26
• @AndrewTawfeek woops, this one slipped under the radar Commented May 2, 2020 at 23:57

The following addresses only initial values only.

If all geodesics are tangent to circle

$$x^2+y^2= r_{min}^2 @ \,r=r_{min}$$

as the Clairaut constant then for an ellipsoid of revolution

$$(z/c)^2+(r/a)^2 = (z/c)^2+(x^2+y^2)/a ^2 = 1$$

we differentiate to find meridian slope

$$\tan \phi=-(a/c) \sqrt{(a/r)^2-1},$$

and $r=f(\theta)$ is found by integration (to get elliptic integral result)

$$\frac{dr}{r \, d\theta}= \sin\phi \, \sqrt{1-(r/r_{min})^2 }.$$

similarly $z$ can be set up. We have $(r,z)$ as functions of $\theta$

With given $(r_B,r_A), \theta_{AB}, (z_B,z_A)$ we cannot directly use this, however..

we can start with initial conditions $(r_B, \theta_A =0, z_A )$ but to arrive at last point we can vary /adjust Clairaut constant $r_{min}$ several small increments in a numerical procedure.

The latest version of geographiclib-octave solves the inverse geodesic problem for a triaxial ellipsoid. The method entails solving the direct problem by simply integrating the ODEs for the geodesic and then using the properties of geodesics on an ellipsoid (as given by Jacobi's solution) to solve the inverse problem. A description of the underlying methods is given here. N.B., for this problem, MATLAB is about 40 times faster than Octave.