Two_body_problem: scipy.integrate.RK45 给出广播错误并且 scipy.integrate.LSODA 从未进入 twoBody 函数

Question

我正在研究二体问题的轨迹计算器，我正在尝试使用 Scipy 的 RK45 or LSODA 来求解 ODE 和 return 轨迹. （如果您认为有 better/easier 方法可以做到这一点，请提出另一种方法）

我将 Spyder IDE 与 Anaconda 一起使用。 Scipy 版本 1.1.0

问题：

RK45：使用 RK45 时，第一步似乎有效。在调试器中单步执行代码时，输入了 twoBody()，并且在第一个运行中完全按照预期工作。然而，在第一个 return ydot 之后，事情就开始出错了。在行 ydot[0] = y[3] 处设置断点，我们开始看到问题所在。数组 y（我预计是 6x1）现在是 6x6 数组。尝试评估此行时，numpy returns 错误 ValueError: could not broadcast input array from shape (6) into shape (1)。我的代码中是否存在导致 y 从 6x1 变为 6x6 的错误？下面是广播错误 returned.

之前的数组 y

y = 
 -5.61494e+06   -2.01406e+06    2.47104e+06     -683.979    571.469  1236.76
-5.61492e+06    -2.01404e+06    2.47106e+06     -663.568    591.88   1257.17
-5.6149e+06     -2.01403e+06    2.47107e+06     -652.751    602.697  1267.99
-5.61492e+06    -2.01405e+06    2.47105e+06     -672.901    582.547  1247.84
-5.61492e+06    -2.01405e+06    2.47105e+06     -672.988    582.46   1247.75
-5.61492e+06    -2.01405e+06    2.47105e+06     -673.096    582.352  1247.64

我的初始条件 Y0 会不会导致它到达的步长太小，因此出错？

苏打水：我还尝试使用 LSODA 求解器。但是，它甚至从未进入 twoBody() 函数！ twoBody() 顶部的断点永远不会到达，程序 return 是运行时间。我不知道这里发生了什么。估计是我设置错了。

编辑：使用 Scipy 的 solve_ivp 时也会发生同样的情况。其他所有方法整合return广播错误

import numpy as np
import scipy.integrate as ode
from time import time
startTime = time()

def twoBody(t, y):
    """
    Two Body function returns the derivative of the state space variables.
INPUTS: 
    --- t ---
        A scalar time value. 

    --- y ---
        A 6x1 array of the state space of a particle in 3D space  
OUTPUTS:
    --- ydot ---
        The derivative of y for the two-body problem

"""
    mu = 3.986004418 * 10**14

    r = np.sqrt(y[0]**2 + y[1]**2 + y[2]**2)

    ydot    = np.empty((6,1))
    ydot[:] = np.nan

    ydot[0] = y[3]             
    ydot[1] = y[4]             
    ydot[2] = y[5]             
    ydot[3] = (-mu/(r**3))*y[0]
    ydot[4] = (-mu/(r**3))*y[1]
    ydot[5] = (-mu/(r**3))*y[2]

    return ydot


# In m and m/s
# first three are the (x, y, z) position
# second three are the velocities in those same directions respectively
Y0 = np.array([-5614924.5443320004,
               -2014046.755686,
               2471050.0114869997,
               -673.03650300000004,
               582.41158099999996,
               1247.7034980000001])

solution = ode.LSODA(twoBody, t0 = 0.0, y0 = Y0, t_bound = 351.0)
#solution = ode.RK45(twoBody, t0 = 0.0, y0 = Y0, t_bound = 351.0)

runTime = round(time() - startTime,6)
print('Program runtime was {} s'.format(runTime))

Answer 1

您也可以使用 scipy.integrate.odeint 来完成这种可能更容易设置的任务：

import numpy as np
from scipy.integrate import odeint


def twoBody(y, t, mu):
    """
    Two Body function returns the derivative of the state space variables.
INPUTS:
    --- t ---
        A scalar time value.

    --- y ---
        A 6x1 array of the state space of a particle in 3D space
OUTPUTS:
    --- ydot ---
        The derivative of y for the two-body problem

"""

    r = np.sqrt(y[0]**2 + y[1]**2 + y[2]**2)

    ydot = np.empty((6,))

    ydot[0] = y[3]
    ydot[1] = y[4]
    ydot[2] = y[5]
    ydot[3] = (-mu/(r**3))*y[0]
    ydot[4] = (-mu/(r**3))*y[1]
    ydot[5] = (-mu/(r**3))*y[2]

    return ydot


# In m and m/s
# first three are the (x, y, z) position
# second three are the velocities in those same directions respectively
Y0 = np.array([-5614924.5443320004,
               -2014046.755686,
               2471050.0114869997,
               -673.03650300000004,
               582.41158099999996,
               1247.7034980000001])

mu = 3.986004418 * 10**14

solution = odeint(twoBody, Y0, np.linspace(0., 351., 100), args=(mu, ))

我无法判断输出是否正确，但似乎集成得很好。

Answer 2

问题是你把ydot定义为矩阵，即二维数组，即使第二维只有宽度1。初始值为一个简单的长度为6的数组。根据Matlab的特点，通过模型自动转换后，y0被解释为行向量，(6,1)数组照常被解释为列向量，求和两者都是 (6,6) 矩阵。

a = np.zeros((3,1))
b = np.array([1,2,3.0])
a+b

array([[1., 2., 3.],
       [1., 2., 3.],
       [1., 2., 3.]])

因此，在下一步中，您将尝试用 (6,6) 矩阵 y 的行填充 (6,1) 列向量 ydot 的条目，这给出报告的错误。

因此避免混合使用不同的矢量类型。定义

ydot = np.empty((6,))

只需稍作改动即可做到这一点。

PS：这可能无关紧要，但我会在 ODE 函数解释之后，直接在求解器调用之前取开始时间。

Two_body_problem: scipy.integrate.RK45 给出广播错误并且 scipy.integrate.LSODA 从未进入 twoBody 函数

Two_body_problem: scipy.integrate.RK45 gives broadcasting error and scipy.integrate.LSODA never enters the twoBody function

python

python-3.x

scipy

ode

runge-kutta