很难为测地方程创建函数 `scipy.integrate.solve_bvp`

Question

自从我开始学习Python库

sympy

和

scipy

以来，我想知道如何使用函数

scipy.integrate.solve_bvp

和

sympy

符号来实现测地方程[1]。我能够为其积分定义实现 Christoffel 符号计算的符号版本。此外，我还定义了积分步骤的 ode 函数。

我无法使模块

scipy.integrate

正常工作，如模块示例 [3] 所示。该库需要一些初始值作为初始猜测。该算法应该收敛到以恒定速度穿过的绑定点的曲线。初始点和终点均未提供，速度范数沿曲线也不是恒定的。

因此请求各位好心人帮忙调试一下这个问题。

[1] https://en.wikipedia.org/wiki/Solving_the_geodesic_equations

[2] https://github.com/alloyha/experiments/blob/main/data/geodesics/ellipsoid_geodesics.ipynb

[3] https://docs.scipy.org/doc/scipy/reference/ generated/scipy.integrate.solve_bvp.html#scipy.integrate.solve_bvp

Answer 1

我会将此问题视为最小化问题而不是微分方程积分。

这是一个示例，我根据路径中等距点的分段线性插值计算椭球上的测地线近似值。

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.sin(theta), b * h * np.cos(theta), c * np.sin(phi)

    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