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

Question

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)

I'm using the Spyder IDE with Anaconda. Scipy version 1.1.0

THE PROBLEMS:

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

Could my initial condition Y0 be causing it to reach a step too small, and therefore error out?

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.

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

I'll take it out and see what happens – natemcintosh Jul 28 '18 at 22:42 — natemcintosh, Jul 28 '18 at 22:42
Having replaced it with np.zeros((6,1)), nothing changes – natemcintosh Jul 28 '18 at 22:45 — natemcintosh, Jul 28 '18 at 22:45

Cleb · Accepted Answer · 2018-07-28T23:29:48.823

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.

score 0 · Answer 2 · answered Feb 03 '19 at 13:48

The problem is that you define ydot as matrix, that is, 2-dimensional array, even if the second dimension only has width 1. The initial value is a simple array of length 6. By automatic conversion after the model of the idiosyncrasies of Matlab, y0 gets interpreted as row vector, the (6,1) array as usual as column vector, the sum of both is a (6,6) matrix.

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

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

So in the next step you are trying to fill the entries of the (6,1) column vector ydot with the rows of the (6,6) matrix y, which gives the reported error.

Thus avoid mixing different vector types. Defining

ydot = np.empty((6,))

does this with a minimal change.

PS: It might not matter, but I would take the start time after the interpretation of the ODE function, directly before the solver call.

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

2 Answers2

Linked