Since I started learning python library sympy and scipy, I wondered how to implement the geodesic equation [1] using the function scipy.integrate.solve_bvp and sympy symbols. I was able to implement a symbolic version of Christoffel's symbols calculation for its integration definition. Also, I defined the ode function for integration steps.
I am not able to make the module scipy.integrate work properly as on the module example [3]. The library requires some initial value as initial guess. The algorithm should converge to the curve which binds the points, which is traversed with a constant velocity. Neither the initial and final points are the provided not the velocity norm is constant along the curve.
Therefore, I ask for the help of some good fellows help debugging this issue.
[1] https://en.wikipedia.org/wiki/Solving_the_geodesic_equations
[2] https://github.com/alloyha/experiments/blob/main/data/geodesics/ellipsoid_geodesics.ipynb
I would approach this problem as a minimization problem not a differential equation integration.
Here is an example where I compute the geodesic approximation on an ellipsoid based on the piece-wise linear interpolation over equally spaced points in the path.
from scipy.optimize import minimize
import numpy
class Ellipsoid:
def __init__(self, a=1, b=1, c=1):
self.abc = [a,b,c]
def path(self, theta, phi):
'''
Path over a paraboloid with phi=0 at the equator,
and at the poles +pi/2 and -pi/2
'''
h = np.cos(phi)
a,b,c = self.abc
return a * h * np.cos(theta), b * h * np.sin(theta), c * np.sin(phi)
def discretized_path_length(self, theta, phi):
'''
Get the points from a path and compute the length
of the piecewise linear interpolation of it
'''
return np.sqrt(sum(np.diff(u)**2 for u in self.path(theta, phi))).sum()
def _discretized_packed_path_length(self, packed_path, p0, p1):
'''
Interface to scipy minimization
take a packed array with the parametized path, and two extra
arguments, p0 the initiala point parameters, and p1 the final
point parameters.
'''
theta, phi = np.vstack([[p0], packed_path.reshape(-1, 2), [p1]]).T
return self.discretized_path_length(theta, phi)
def geodesic(self, p1, p2, n, method='least_squares'):
'''
Given p1=(theta1, phi1), p2=(theta1, phi2) and a number
of segments computes the discretized geodesic
'''
theta1, phi1 = p1
theta2, phi2 = p2
t_theta = np.linspace(theta1, theta2, n+1)
t_phi = np.linspace(phi1, phi2, n+1)
t_packed = np.array([t_theta, t_phi]).T
results = minimize(
fun=self._discretized_packed_path_length,
x0=t_packed[1:-1].reshape(-1),
args=(t_packed[0], t_packed[-1]))
t_packed[1:-1] = results.x.reshape(-1, 2)
theta, phi = t_packed.T
return theta, phi
To validate the approach we can check the case where the ellipsoid reduces to a sphere and the geodesic length can be easily calculated as the length of the arc connecting two points.
def check_geodesic(p0, p1, n):
ball = Ellipsoid(1,1,1)
theta, phi = zip(p0, p1)
theta0 = np.linspace(*theta, n+1)
phi0 = np.linspace(*phi, n+1)
initial_length = ball.discretized_path_length(theta0, phi0)
theta, phi = ball.geodesic(p0, p1, n, 'minimize')
m_geodesic_length = ball.discretized_path_length(theta, phi)
distance = ball.discretized_path_length(*zip(p0, p1))
arc = 2*np.arcsin(distance/2)
print(theta[0], phi[0], theta[-1], phi[-1])
print('Initial path length:', initial_length)
print('Distance from start to end', distance)
print('Geodesic length (minimize):', m_geodesic_length)
print('Length of arc from start to end', arc)
A trivial case traveling over the equator
check_geodesic((0, 0), (1, 0), 100)
0.0 0.0 1.0 0.0
Initial path length: 0.9999958333385416
Distance from start to end 0.9588510772084059
Geodesic length (minimize): 0.9999958333385416
Length of arc from start to end 0.9999999999999999
A non-trivial geodesic
check_geodesic((0, 0.5), (1, 0.5), 100)
0.0 0.5 1.0 0.5
Initial path length: 0.8775789053009355
Distance from start to end 0.8414709848078966
Geodesic length (minimize): 0.8685090778499669
Length of arc from start to end 0.8685118212476727