代码之家 › 专栏 › 技术社区 › natemcintosh

二体问题:scipy.integrate.rk45给出广播错误,scipy.integrate.lsoda从不进入two body函数

runge-kutta ode scipy python-3.x python

natemcintosh · 技术社区 · 6 年前

我正在做一个双体问题的轨迹计算器,我正在尝试使用scipy的 RK45 或 LSODA 解出ODE并返回轨迹。(如果您认为有更好/更简单的方法,请建议其他方法)

我在用刺探器和水蟒。scipy版本1.1.0

问题:

RK45: 使用时 RK45 第一步似乎可行。在调试程序中单步执行代码时,会输入twobody(),并在第一次运行时完全按预期工作。但是,在第一次 return ydot ,事情开始出错。行上有断点 ydot[0] = y[3] 我们开始发现问题所在。数组 y (我希望是6x1)现在是6x6数组。尝试计算此行时,numpy返回错误 ValueError: could not broadcast input array from shape (6) into shape (1) . 我的代码中是否有可能导致 Y 从6x1到6x6?下面是返回广播错误前的数组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() 功能!顶部的内部断点 二体() 从未到达,程序返回运行时。我不知道这是怎么回事。我猜我设置不正确。

编辑: 同样的情况也发生在使用scipy's时。 solve_ivp . 所有其他的集成方法都返回广播错误。

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

2 回复 | 直到 5 年前

Cleb 6 年前

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

Lutz Lehmann 5 年前

ydot 1 y0

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

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

y

ydot = np.empty((6,))