如何用Python脚本数值求解二阶二元PDE边值问题并修正代码？

问题背景 #

我有如下两个二阶微分方程（伪代码）：

r'' = (s')² * r - G*M*r⁻²
r' = -(r*s'')/(2*s')

其中s和r是关于时间的函数，'表示对时间的导数，G和M为常数。这组方程来自月球绕地球的二维极坐标Euler-Lagrange方程，方程正确性无需质疑。

我尝试用scipy求解，但该库要求先把方程拆分为仅含r或s的一阶ODE，我不知道怎么操作，而且复杂版本的方程拆分难度很大。想问：

• 有没有办法无需任何变换直接求解这类方程？
• 是否有相关库支持这种求解方式？

另外这是边值问题：已知时间区间首尾的r和s值（来自JPL星历），但未知r'和s'的初值。

代码问题与修正 #

参考建议后我写了一段scipy代码，但无论时间区间长短，最终theta都会变成0，显然结果不对，而且方程在equations函数里的写法肯定有误，请问怎么正确表示？

错误代码 #

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from skyfield.api import load, Topos

G = 6.67430e-11
M = 5.97219e24

eph = load("de421.bsp")
observer = eph["earth"] + Topos(latitude_degrees=0, longitude_degrees=0)
moon = eph["moon"]

time1 = load.timescale().utc(2000, 1, 1, 0, 0, 0)
moon_position = observer.at(time1).observe(moon)
theta0, _, r0 = moon_position.radec()
time2 = load.timescale().utc(2000, 1, 1, 0, 0, 1)
moon_position = observer.at(time2).observe(moon)
theta1, _, r1 = moon_position.radec()

theta0 = float(theta0.radians)
thetadot = theta0 - float(theta1.radians)
r0 = float(r0.m)
rdot = r0 - float(r1.m)
print(r0, r1.m)
print(theta0, theta1.radians)
print(thetadot, rdot)


def equations(t, y):
    theta, r, thetadot, rdot = y
    dydt = [rdot, (thetadot**2 * r - G * M / r**2), thetadot, -r * (thetadot**2) / (2 * rdot),]
    return dydt


t_eval = np.linspace(0, 10000, 1000000)
sol = solve_ivp(equations, [0, 10000], [r0, theta0, rdot, thetadot])
t_plot = sol.t
y_plot = sol.y

#plt.plot(t_plot, y_plot[0], label="r(t)")
plt.plot(t_plot, y_plot[1], label="theta(t)")
plt.xlabel("Time")
plt.ylabel("Values")
plt.legend()
plt.show()

错误分析与修正方案 #

1. 状态变量与导数对应完全错误

你的状态变量y定义为theta, r, thetadot, rdot，但导数列表dydt的顺序完全不符合逻辑：

• theta的导数应为thetadot，而非rdot
• r的导数应为rdot，而非第一个方程的r''
• thetadot的导数是s''（代码中用theta替代了原方程的s），需要从第二个原方程推导
• rdot的导数是r''，对应第一个原方程

2. 原方程的正确转换

将原方程中的s替换为代码里的theta，整理后：

• 方程1：r'' = (theta')² * r - G*M/r² → 这是rdot的导数
• 方程2：r' = -(r*theta'')/(2*theta') → 变形得theta'' = -(2*r'*theta')/r，这是thetadot的导数

3. 修正后的equations函数

def equations(t, y):
    theta, r, thetadot, rdot = y
    # 基础导数
    dtheta_dt = thetadot
    dr_dt = rdot
    # 推导theta''（避免除零错误）
    dthetadot_dt = -(2 * rdot * thetadot) / r if r != 0 else 0
    # 推导r''（避免除零错误）
    drdot_dt = (thetadot**2 * r) - (G * M) / (r**2) if r != 0 else 0
    return [dtheta_dt, dr_dt, dthetadot_dt, drdot_dt]

4. 初值计算错误

你当前计算thetadot和rdot时，用了反向差商且未除以时间间隔，导致角速度为负，最终theta不断减小到0。正确的差商计算应基于时间间隔（这里为1秒）：

delta_t = 1.0  # time1到time2的时间间隔
thetadot = (float(theta1.radians) - theta0) / delta_t  # 正向差商，符合月球赤经递增规律
rdot = (float(r1.m) - r0) / delta_t

5. 边值问题的正确求解方式

solve_ivp是用于初值问题的求解器，你的问题是边值问题，需要用scipy.integrate.solve_bvp，示例框架如下：

from scipy.integrate import solve_bvp

def bvp_equations(t, y):
    theta, r, thetadot, rdot = y
    dtheta_dt = thetadot
    dr_dt = rdot
    dthetadot_dt = -(2 * rdot * thetadot) / r
    drdot_dt = (thetadot**2 * r) - (G * M) / (r**2)
    return np.vstack([dtheta_dt, dr_dt, dthetadot_dt, drdot_dt])

def bvp_boundary(ya, yb):
    # 边界条件：匹配首尾的theta和r值
    return np.array([ya[0] - theta0, yb[0] - theta1, ya[1] - r0, yb[1] - r1])

# 初始化状态猜测
t_guess = np.linspace(0, 10000, 100)
y_guess = np.zeros((4, t_guess.size))
y_guess[0] = np.linspace(theta0, theta1, t_guess.size)
y_guess[1] = np.linspace(r0, r1, t_guess.size)
y_guess[2] = thetadot  # 初始角速度猜测值
y_guess[3] = rdot      # 初始径向速度猜测值

sol_bvp = solve_bvp(bvp_equations, bvp_boundary, t_guess, y_guess)

关于无需变换直接求解的问题 #

目前主流Python ODE求解库（如scipy）几乎都要求将高阶ODE转换为一阶方程组，这是数值求解的标准流程——所有高阶ODE的数值算法本质上都是基于一阶方程组设计的。

如果不想手动拆分，可以用sympy等符号计算库自动生成一阶方程组：

import sympy as sp

t = sp.symbols('t')
r = sp.Function('r')(t)
s = sp.Function('s')(t)
G, M = sp.symbols('G M')

# 定义原方程
eq1 = sp.Eq(r.diff(t,2), (s.diff(t))**2 * r - G*M / r**2)
eq2 = sp.Eq(r.diff(t), -(r * s.diff(t,2))/(2*s.diff(t)))

# 定义一阶变量
r1 = r.diff(t)
s1 = s.diff(t)
# 自动求解二阶导数
sol = sp.solve([eq1, eq2], [r.diff(t,2), s.diff(t,2)])
# 得到一阶方程组：r'=r1, r1'=sol[r''], s'=s1, s1'=sol[s'']

之后可将符号表达式转换为数值函数，再用scipy求解。

内容的提问来源于stack exchange，提问作者AtomProgrammer