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
问题描述
我正在为两体问题"开发轨迹计算器，并且尝试使用Scipy的 LSODA 来求解ODE并返回轨迹. (如果您认为这样做的更好/更简便的方法，请提出另一种方法)
I'm working on a trajectory calculator for the Two Body Problem, and I'm attempting to use Scipy's RK45 or LSODA to solve the ODE and return the trajectory. (Please suggest another method if you think there's a better/easier way to do this)
我正在将Spyder IDE与Anaconda一起使用. Scipy版本1.1.0
I'm using the Spyder IDE with Anaconda. Scipy version 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.
RK45:
When using RK45, the first step seems to work. When stepping through the code in the debugger, twoBody() is entered, and works exactly as expected the first run through. However, after the first return ydot, things start to go wrong. With a breakpoint at the line ydot[0] = y[3], we start to see the problem. The array y (which I expected to be 6x1) is now a 6x6 array. When attempting to evaluate this line, numpy returns the error
ValueError: could not broadcast input array from shape (6) into shape (1). Is there an error in my code that would cause y to go from 6x1 to 6x6? Below is the array y right before the broadcasting error is returned.
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会导致它到达的步幅太小，从而导致错误吗?
Could my initial condition Y0 be causing it to reach a step too small, and therefore error out?
LSODA:
我还尝试使用 LSODA 求解器.但是，它甚至都不会进入twoBody()函数！永远不会到达twoBody()顶部的内部断点，并且程序将返回运行时.我不知道这是怎么回事.猜猜我设置不正确.
LSODA:
I also tried to use the LSODA solver. However, it never even enters the twoBody() function! A breakpoint inside at the top of twoBody() is never reached, and the program returns the runtime. I have no idea what's going on here. Guessing I set it up incorrectly.
使用Scipy的 resolve_ivp .所有其他集成方法都将返回广播错误.
The same occurs when using Scipy's solve_ivp. All the other methods of integration return the broadcast error.
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))
推荐答案
You can also use scipy.integrate.odeint for this kind of task which might be easier to set up:
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, ))
我无法判断输出是否正确，但似乎集成得很好.
I cannot judge whether the output is correct, but it seems to integrate well.
这篇关于Two_body_problem:scipy.integrate.RK45给出广播错误，并且scipy.integrate.LSODA从未进入twoBody函数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
