python is freezing when executing a 3d plot

Question

From plotting orbital trajectories, we have the following code. The values have been changed to known IC that work.

If this code is correct (it can't be though), it would generate

Running this code just freeze my computer or outputs absolute wrong plots. Could someone help me find a way to fix this?

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from numpy import linspace
from mpl_toolkits.mplot3d import Axes3D

mu = 398600
# r0 = [-149.6 * 10 ** 6, 0.0, 0.0]  #  Initial position
# v0 = [29.9652, -5.04769, 0.0]      #  Initial velocity
u0 = [-4069.503, 2861.786, 4483.608, -5.114, -5.691, -5]


def deriv(u, dt):
    n = -mu / np.sqrt(u[0] ** 2 + u[1] ** 2 + u[2] ** 2)
    return [u[3],     #  dotu[0] = u[3]'
            u[4],     #  dotu[1] = u[4]'
            u[5],     #  dotu[2] = u[5]'
            u[0] * n,       #  dotu[3] = u[0] * n
            u[1] * n,       #  dotu[4] = u[1] * n
            u[2] * n]       #  dotu[5] = u[2] * n

dt = np.arange(0.0, 24 * 3600, .01)   #  Time to run code in seconds'
u = odeint(deriv, u0, dt)
x, y, z, x2, y2, z2 = z.T

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(x, y, z)
plt.show()

Answer 1

You are running into numerical problems for several reasons: Firstly, you should not ask the ODE solver to return data for 8640000 points. Secondly, your parameters and initial conditions contain large numbers, which you could probably get rid of by introducing appropriate non-dimensional quantities.

That being set, the code below produces sensible output:

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from numpy import linspace
from mpl_toolkits.mplot3d import Axes3D

u0 = [0, 0, 1, 0, -1, 0]
mu = .1


def deriv(u, dt):
    n = -mu / np.sqrt(u[0] ** 2 + u[1] ** 2 + u[2] ** 2)
    return [u[3],     #  dotu[0] = u[3]'
            u[4],     #  dotu[1] = u[4]'
            u[5],     #  dotu[2] = u[5]'
            u[0] * n,       #  dotu[3] = u[0] * n
            u[1] * n,       #  dotu[4] = u[1] * n
            u[2] * n]       #  dotu[5] = u[2] * n

times = np.linspace(0.0, 200, 100)
u = odeint(deriv, u0, times)
x, y, z, x2, y2, z2 = u.T

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(x, y, z)
plt.show()

The result is