我会把这个问题作为一个最小化问题来处理，而不是微分方程的积分。
这里有一个例子，我在一个椭圆体上计算测地线近似，基于路径中等距点上的分段线性插值。

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

字符串
为了验证这种方法，我们可以检查椭球体缩小为球体的情况，测地线长度可以很容易地计算为连接两点的弧的长度。

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)

型
一个在赤道上旅行的小例子

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

的字符串
非平凡测地线

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

的字符串