如何正确使用 scipy.integrate 中的 ode() 处理二阶方程？

Question

这些是我的星球的运动方程。 使用odeint，它工作得很好，因为odeint可以处理二阶。 但 odeint 不是很精确，所以我被要求使用 ode 代替。 但是 ode 无法处理二阶方程，您需要将它们转换为一阶方程。 但在这种情况下似乎很困难，因为 d²x/dt² 和 d²ydt² 与 dx/dt 和 dy/dt 无关。

https://paste.pythondiscord.com/EFG5W6VL7JWDHU4TVA4USWFECY 这是我的工作代码，使用 odeint，以及我的非工作代码，使用 ode。 （在我们的例子中仅导入了前 70 行）。

我尝试将 d²x/dt² 和 d²ydt² 转换为 dx/dt 和 dy/dt，将它们乘以 t，然后加上 k1 和 k2（初始速度），但它不起作用。

有什么想法吗？

感谢您的宝贵时间。

Answer 1

但是 ode 无法处理二阶方程

ode

和

odeint

都不是设计用于二阶方程组的。他们都期望方程能够以一阶形式表达。

你需要将这些转换为一阶。

不，你不知道。

您链接到的方程是：

这些已经是一阶形式了。当然，它们代表两个二阶运动方程的系统。例如，您可以将前两个方程组合起来得到二阶方程：

但正如所写，它是一个由四个一阶运动方程组成的系统。

如果

odeint

的工作代码是：

import matplotlib.pyplot as plt
from scipy import integrate
def dzdt(z, t):
    x, vx, y, vy = z
    dxdt = vx
    dvxdt = -x / (x**2 + y**2)**(3/2)
    dydt = vy
    dvydt = -y / (x**2 + y**2)**(3/2)
    return [dxdt, dvxdt, dydt, dvydt]

x0 = 1
vx0 = 0
y0 = 0
vy0 = 1
z0 = [x0, vx0, y0, vy0]
t = np.linspace(0, 10)
res = integrate.odeint(dzdt, z0, t)
plt.plot(res[:, 0], res[:, 1])
plt.axis('square')

您可以获得

ode

的工作代码，例如：

def dzdt2(t, z):
    # reverse the order of the inputs
    return dzdt(z, t)

ode = integrate.ode(dzdt2)
ode.set_integrator('dop853')
ode.set_initial_value(z0, 0)

z = np.zeros((len(t), 4))
for i, ti in enumerate(t):
    z[i] = ode.integrate(ti).real
plt.plot(z[:, 0], z[:, 1])
plt.axis('square')