用 scipy odeint 在 python 中求解向量常微分方程

Posted 2023-03-12

【中文标题】用 scipy odeint 在 python 中求解向量常微分方程

【英文标题】：Solving vector ordinary differential equations in python with scipy odeint

【发布时间】：2015-12-13 20:16:21

【问题描述】：

我正在尝试解决和颂歌涉及向量的问题，但无法得出可行的答案。所以我把它分成6个分量，一个分量的每个时间导数，一个速度分量的每个时间导数。第一个值似乎是合理的，但随后它会跳到数百万的数字，我不知道为什么。老实说，我真的不确定如何首先做到这一点，现在只是试一试。我似乎无法在网上找到任何有关它的信息，如果有此类问题的示例，我可以使用一些帮助或一些链接。任何有关如何解决 ODE 的信息都将不胜感激。

def dr_dt(y, t):
"""Integration of the governing vector differential equation.
d2r_dt2 = -(mu/R^3)*r with d2r_dt2 and r as vecotrs.
Initial position and velocity are given.
y[0:2] = position components
y[3:] = velocity components"""

    G = 6.672*(10**-11)
    M = 5.972*(10**24)
    mu = G*M
    r = np.sqrt(y[0]**2 + y[1]**2 + y[2]**2) 

    dy0 = y[3]
    dy1 = y[4]
    dy2 = y[5]
    dy3 = -(mu / (r**3)) * y[0]
    dy4 = -(mu / (r**3)) * y[1]
    dy5 = -(mu / (r**3)) * y[2]
    return [dy0, dy3, dy1, dy4, dy2, dy5]

解决这个问题后，我想绘制它。它应该变成一个椭圆，但老实说，我也不完全确定如何做到这一点。我正在考虑获取位置的大小，然后随时间绘制它。如果有更好的方法，请随时告诉我。

谢谢。

【参考方案1】：

首先，如果您的y[:3] 是位置并且y[3:] 是速度，那么dr_dt 函数应该完全按照这个顺序返回组件。其次，要绘制轨迹，我们可以使用出色的 matplotlib mplot3d 模块，或者省略位置和速度的 zth 分量（因此我们的运动在 XY 平面上），并绘制 y 与 x。示例代码（带有更正的返回值顺序）如下：

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def dr_dt(y, t):
    """Integration of the governing vector differential equation.
    d2r_dt2 = -(mu/R^3)*r with d2r_dt2 and r as vecotrs.
    Initial position and velocity are given.
    y[0:2] = position components
    y[3:] = velocity components"""

    G = 6.672*(10**-11)
    M = 5.972*(10**24)
    mu = G*M
    r = np.sqrt(y[0]**2 + y[1]**2 + y[2]**2).

    dy0 = y[3]
    dy1 = y[4]
    dy2 = y[5]
    dy3 = -(mu / (r**3)) * y[0]
    dy4 = -(mu / (r**3)) * y[1]
    dy5 = -(mu / (r**3)) * y[2]
    return [dy0, dy1, dy2, dy3, dy4, dy5]

t = np.arange(0, 100000, 0.1)
y0 = [7.e6, 0., 0., 0., 1.e3, 0.]
y = odeint(dr_dt, y0, t)
plt.plot(y[:,0], y[:,1])
plt.show()

这会产生一个漂亮的椭圆：

【讨论】：

非常感谢，这很有帮助。我想我可能会把这些弄乱，但主要问题是当我注意到你将y0 保留为列表时。我有我的数组，当我发布它时，我应该把它放在代码中，我的错误。一旦我修复它看起来好像更多相关信息正在返回。我现在唯一的问题是它无法正确绘制。我有y0 = [-4069.503, 2861.786, 4483.608, -5.114, -5.691, -1.500] 供输入。当它返回时，我绘制它。它只是从左上角到右下角的一条线。当我尝试时，3D 也是如此。

由于您的计算采用国际单位制，因此您的初始速度需要以 m/s 为单位（以及分别以 m 为单位的位置）。 5 m/s 的速度远不及逃逸速度，因此您可以看到物体落到地球上的轨迹。

难以置信，完全有效。非常感谢，我不敢相信这是一个如此简单的错误。非常感谢您的帮助。