如何使用 python 脚本对一组两个二阶双变量偏微分方程进行数值近似解（边界值问题）？

Question

我有一组两个偏微分方程，如下所示（伪代码）：

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

s

和

r

是时间的函数，

'

是导数，

G

和

M

是常数。

它们是从二维极坐标中月球绕地球运动的欧拉-拉格朗日方程获得的。 他们甚至可能是错的；他们的正确性不是这个问题的重点。

我尝试使用

scipy

库来求解它们，但它要求我首先将方程分成一阶 ODE，以便每个方程包含 r 或 s，但不能同时包含两者。坦率地说，我不知道该怎么做，虽然这对于这两个方程来说是可能的，但对于我计划使用的更复杂的版本来说，它不太可能实现。

所以我的问题是：我怎样才能在不进行任何操作的情况下直接解决它们？有没有图书馆可以做到这一点？

此外，这是一个边值问题。我知道给定时间间隔开始和结束时的 r 和 s 值（来自 JPL 星历），但不知道

r'

或

s'

。

编辑： 看了评论后，我用 scipy 组合了一些代码：

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()

不幸的是，无论我的时间尺度多长或多短，theta 总是在最后变成 0，这是不正确的。另外，我在方程函数中表示方程组的方式是错误的，我该如何正确地做到这一点？

Answer 1

为了尝试找到解析解（尽管我注意到该问题要求数值近似），您可以尝试使用 sympy，特别是

dsolve_system

函数，例如，

from sympy import *
from sympy.solvers.ode.systems import dsolve_system

s, r, G, M, t = symbols("s r G M t")

sf = Function(s)(t)
rf = Function(r)(t)

sp = Derivative(sf, t)  # derivative
spp = Derivative(sf, (t, 2))  # second derivative

rp = -(r * spp) / (2 * sp)
rpp = sp**2 * r - G * M / r**2

eqs = [Eq(Derivative(rf, t), rp), Eq(Derivative(rf, (t, 2)), rpp)]

dsolve_system(eqa, t=t)

不幸的是，sympy 似乎无法解析求解这个方程组。