[Tutor] Problems using odeint

[eryksun]

On Tue, Jan 14, 2014 at 1:40 PM, Grace Roberts
<xgrace_robertsx at hotmail.co.uk> wrote:
>
> -For certain initial conditions the programme displays impossible orbits,
> showing the satellite making immediate sharp turns or jumping up and down in
> velocity. The exact same code can also end up displaying different graphs
> when run multiple times. For example when initial conditions are set at:
> xx0=[1000.,10000.,1000.,10000.].

The trajectory is decaying into the planet. In real life it hits the
surface. In the model it continues until the radius is 0. As the
radius approaches zero, odeint will have problems converging.

I ran your calculation for Phobos, assuming a perfectly circular
orbit. I started with `Rsx == Rxy` and `Vsx = -Vsy`. It steps through
a cycle in the expected time.

https://en.wikipedia.org/wiki/Phobos_%28moon%29

Output Plot:
http://i.imgur.com/Jh0sbnT.png


    import numpy as np
    from scipy import integrate
    import matplotlib.pyplot as plt

    G = 6.674e-11
    Mm = 6.4185e23  # Mars mass
    Rs = 9.376e6    # Phobos orbit semi-major axis
    Vs = 2.138e3    # Phobos orbit average velocity
    Ts = 2.7554e4   # Phobos orbit period

    def f(y, t):
        rx, ry, vx, vy = y
        r = np.hypot(rx, ry)
        ax = -G * Mm * rx / r ** 3
        ay = -G * Mm * ry / r ** 3
        return vx, vy, ax, ay

    Rsx = Rsy = (Rs ** 2 / 2) ** 0.5
    Vsy = (Vs ** 2 / 2) ** 0.5
    Vsx = -Vsy
    y0 = [Rsx, Rsy, Vsx, Vsy]

    t = np.linspace(0, Ts, Ts * 10)

    y, info = integrate.odeint(f, y0, t, full_output=1)
    rx, ry, vx, vy = y.T

    plt.subplot(211)
    plt.title('Position')
    plt.plot(t, rx, t, ry)
    plt.grid(); plt.axis('tight')
    plt.subplot(212)
    plt.title('Velocity')
    plt.plot(t, vx, t, vy)
    plt.grid(); plt.axis('tight')
    plt.show()