Solving vector ordinary differential equations in python with scipy odeint
Question: I'm trying to solve and ode that involves vectors and can't come out with a
feasible answer. So i split it into 6 components, one for each time derivative
of a component, and one for each time derivative of a velocity component. The
first value seems reasonable but then it jumps to to numbers in the millions
and I'm not sure why. I'm honestly not really sure how to do this in the first
place and am just giving it a shot right now. I couldn't seem to find any info
on it online and could use some help or some links if there are examples of
this type of problem. Any info would be much appreciated on how to get this to
solve the 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]
    

After this is solved, I want to plot it. It should come out to an ellipse but
to be honest I'm not exactly sure how to do that either. I was thinking of
taking the magnitude of the position and then plotting it with time. If
there's a better way of doing this please feel free to let me know.

Thanks.


Answer: First, if your `y[:3]` is position and `y[3:]` is velocity, then `dr_dt`
function should return the components in exactly this order. Second, to plot
the trajectory we can either use excellent matplotlib `mplot3d` module, or
omit `z`th component of position and velocity (so our motion is on XY plane),
and plot `y` versus `x`. The sample code (with corrected order of return
values) is given below:

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

This yields a nice ellipse:

[![enter image description
here](http://i.stack.imgur.com/tey09.png)](http://i.stack.imgur.com/tey09.png)

