Jacobi's solution for geodesics on a triaxial ellipsoid

✨

This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine the Earth not as a perfect sphere or even a slightly squashed ball (like a standard basketball), but as a lumpy, three-sided potato. In math terms, this is a triaxial ellipsoid. It has three different lengths: a long side, a medium side, and a short side.

For centuries, scientists knew how to draw the "straightest possible line" (a geodesic) between two points on a sphere or a standard squashed ball. But on this lumpy, three-sided potato, the math was a nightmare. It was like trying to navigate a maze where the walls keep shifting shape.

In 1838, a brilliant mathematician named Carl Gustav Jacob Jacobi found a secret key. He discovered that even on this lumpy potato, you could describe the path using simple one-dimensional integrals (basically, summing up tiny slices of a curve). However, for nearly 200 years, nobody knew how to actually use his formula on a computer without it crashing or giving wrong answers.

Charles Karney, the author of this paper, finally cracked the code. Here is how he did it, explained simply:

1. The Problem: The "Lumpy Potato" Navigation

If you try to walk in a straight line on a lumpy potato, your path gets weird. It might loop around the long way, or zigzag back and forth.

The Old Way: Previous methods tried to simulate the walk step-by-step (like a video game character taking tiny steps). This works for short trips, but if you walk for a long time, the tiny errors pile up, and you end up in the wrong country.
The Goal: Karney wanted a method that is accurate no matter how far you walk, and fast enough to be useful.

2. The Solution: The "Fourier Fourier" Trick

Jacobi's solution involves complex math formulas (integrals) that are hard to calculate directly. Karney's big idea was to treat these formulas like a song.

The Analogy: Imagine the complex curve of the path is a complicated musical note. Instead of trying to calculate the exact physics of the sound wave, Karney broke the note down into a Chorus of Sine Waves (a Fourier series).
Why it works: Sine waves are easy to integrate (sum up). By approximating the difficult math as a sum of simple waves, he could calculate the path with extreme precision. It's like turning a jagged, rocky mountain path into a smooth, rolling hill that a computer can measure instantly.

3. The "Umbilic" Trap

There are two special spots on the potato (called umbilics) where the surface is perfectly round, like the top of a sphere.

The Problem: At these spots, the math formulas break down and scream "Error!" (mathematical singularities).
The Fix: Karney realized that geodesics passing through these spots behave in a predictable, chaotic way. He created a special "bridge" to jump over these spots, calculating exactly how the path twists and turns as it passes through the roundest part of the potato.

4. Solving the "Inverse" Puzzle

There are two types of navigation problems:

Direct: "I am here, facing North. Where do I end up after 100 miles?" (Easy with Karney's method).
Inverse: "I am here, and I want to go to that specific point. Which direction should I face?" (Harder).

Karney used a clever trick for the second problem. He treated the direction you face as a dial. He spun the dial, calculated where you would end up, and kept adjusting the dial until you landed exactly on your target. Because his "song" method was so accurate, this dial-spinning process converged (found the answer) very quickly.

5. Why This Matters

For Earth: While Earth is mostly a standard squashed ball, some parts of it (and other planets like Mars or asteroids) are lumpy. This method allows for ultra-precise mapping of these weird shapes.
For Space: If you are sending a probe to an asteroid that looks like a potato, you need to know exactly how to fly around it without crashing.
For Math: It proves that Jacobi's 180-year-old secret was correct and usable, turning a theoretical curiosity into a practical tool.

The Bottom Line

Charles Karney took a 19th-century mathematical masterpiece that was too messy to use, cleaned it up with modern digital tools (Fourier series), and built a high-precision GPS for lumpy, three-sided worlds. He showed that even on the most irregular shapes, the "straightest path" can be found with incredible speed and accuracy.

1. Problem Statement

The paper addresses the geodesic problem on a triaxial ellipsoid (an ellipsoid with three distinct semi-axes $a > b > c$ ). While geodesics on biaxial ellipsoids (revolution ellipsoids, modeling the Earth) are well-understood and solvable via elliptic integrals due to rotational symmetry (Clairaut's constant), the triaxial case lacks obvious symmetry.

Direct Problem: Given a starting point, an initial azimuth, and a distance, find the endpoint and final azimuth.
Inverse Problem: Given two points, find the shortest path (distance and initial/final azimuths) connecting them.
Challenge: Previous efforts to implement Jacobi's 1838 analytical solution (which reduces the problem to one-dimensional integrals) were lacking. Numerical integration of Ordinary Differential Equations (ODEs) is an alternative but suffers from error accumulation over long distances and difficulty in solving the inverse problem.

2. Methodology

Karney implements a numerical solution based on Jacobi's 1839 solution, which utilizes ellipsoidal coordinates ( $\beta, \omega$ ) to separate variables. The core methodology involves:

A. Mathematical Formulation

Ellipsoidal Coordinates: The surface is parameterized by latitude $\beta$ and longitude $\omega$ . The metric is expressed in a form satisfying the "Liouville surface" condition, allowing the geodesic equations to be decoupled.
Conserved Quantity: A constant $\gamma$ $γ$ (generalizing the Clairaut constant) determines the geodesic's behavior.
- $\gamma > 0$ : Circumpolar geodesics (oscillate in $\beta$ , rotate in $\omega$ ).
- $\gamma < 0$ : Transpolar geodesics (oscillate in $\omega$ , rotate in $\beta$ ).
- $\gamma = 0$ : Umbilical geodesics (pass through the umbilic points where principal curvatures are equal).
Integral Representation: The path and distance are given by coupled equations involving integrals of the form:
$\delta = \int f(\beta) d\beta \mp \int f(\omega) d\omega$
$s = \int g(\beta) d\beta \pm \int g(\omega) d\omega$

B. Numerical Implementation Strategy

To achieve high accuracy and efficiency, the paper introduces several specific numerical techniques:

Fourier Series Approximation: Instead of direct numerical quadrature, the integrands are approximated by Fourier series. The coefficients are computed using the Fast Fourier Transform (FFT).
- This allows for the rapid evaluation of indefinite integrals (essential for finding arbitrary points along a geodesic).
- The method supports arbitrary precision; the number of Fourier coefficients scales linearly with the required precision bits.
Variable Transformations: To handle near-singularities (e.g., when $\gamma \to 0$ or near umbilics), the integration variables are transformed using Jacobi elliptic functions (specifically the amplitude function am). This smooths the integrands, preventing the need for excessive Fourier terms.
Root Finding for the Inverse Problem:
- The inverse problem is reduced to a 1D root-finding problem: finding the initial azimuth $\alpha_1$ such that the geodesic hits the target longitude $\omega_2$ for a given latitude $\beta_2$ .
- The author uses Chandrupatla's method (a robust hybrid of bisection and inverse quadratic interpolation) to solve this.
2D Newton's Method with Bisection: For the direct problem (solving for endpoint coordinates given distance), a 2D Newton's method is employed. To ensure convergence, the author implements a bounding box strategy (bisection) that updates based on the signs of the residual functions, preventing the solver from diverging near singularities.

C. Special Cases

Umbilical Geodesics ( $\gamma=0$ ): These are treated as a limiting case. The solution involves tracking the geodesic through multiple passages of the umbilics using a connection relation that relates the azimuth before and after crossing an umbilic.
Biaxial Limit: The implementation naturally reduces to the standard biaxial solution when $a=b$ or $b=c$ , recovering known results.

3. Key Contributions

First Robust Implementation: This is the first comprehensive, high-accuracy numerical implementation of Jacobi's solution for the triaxial ellipsoid.
High-Precision Arithmetic: The algorithm is designed to work with arbitrary precision (beyond standard double precision), allowing users to trade CPU time for accuracy.
Unified Framework: It provides a single framework that handles the direct and inverse problems, as well as the special cases of umbilical geodesics and the biaxial limit.
GeographicLib Integration: The code is released as part of GeographicLib (v2.7), making these advanced geodesic calculations accessible to the scientific community.

4. Results and Performance

The author tested the implementation on "Cayley's ellipsoid" ( $a=\sqrt{2}, b=1, c=1/\sqrt{2}$ ) and other shapes, including a triaxial Earth model.

Accuracy:
- Direct Problem: Mean error of 5 ulp (units in the last place) for position and 6 ulp for direction.
- Inverse Problem: Mean error in distance of 3 ulp.
- Maximum errors observed were roughly 160 ulp (position) and 1500 ulp (direction), which remain small but highlight the difficulty of points near umbilics.
Convergence:
- The 2D Newton method for the direct problem exhibits quadratic convergence.
- The inverse problem solver converges with a rate where the log-error is multiplied by $\approx 0.75$ per iteration (slower than quadratic but robust).
Performance:
- Direct Problem: ~53 $\mu$ s per solution (on an Intel i7).
- Inverse Problem: ~220 $\mu$ s per solution.
- Comparison: The triaxial solution is roughly 10 times slower than the biaxial solution but maintains high accuracy.
Comparison with ODE Integration:
- While ODE integration (e.g., Runge-Kutta) is faster for short distances, its error grows linearly with distance.
- Jacobi's solution maintains constant accuracy regardless of distance, making it superior for long geodesics.

5. Significance

Geodesy and Planetary Science: This work enables precise modeling of shortest paths on bodies that are not rotationally symmetric, such as certain asteroids, moons, or a more accurate model of the Earth (which is slightly triaxial).
Mathematical Physics: It demonstrates the practical utility of 19th-century analytical mechanics (Jacobi's separation of variables) when combined with modern numerical analysis (Fourier series and high-precision arithmetic).
Algorithmic Robustness: The paper provides a blueprint for solving complex, coupled non-linear systems involving singular integrals, offering a stable alternative to direct ODE integration for global pathfinding problems.

In summary, Karney successfully bridges the gap between Jacobi's theoretical solution and practical, high-precision computation, providing a reliable tool for geodesy on triaxial ellipsoids.