Building orbital trajectories in python

Question

How to configure three body problems in python? How to define a function to solve ODE?

Three equations: x'' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * x ,
y'' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * y and
z'' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * z .

Written as 6 first order, we have

x' = x2 ,

y' = y2 ,

z' = z2 ,

x2' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * x ,

y2' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * y and

z2' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * z

I also want to add to the Plot o path the orbit of the Earth and Mars, which we can consider circular. The Earth is located at a distance of 149.6 * 10 ** 6 km from the Sun and Mars 227.9 * 10 ** 6 km.

 #!/usr/bin/env python # This program solves the 3 Body Problem numerically and plots the trajectories import pylab import numpy as np import scipy.integrate as integrate import matplotlib.pyplot as plt from numpy import linspace mu = 132712000000 #gravitational parameter r0 = [-149.6 * 10 ** 6, 0.0, 0.0] v0 = [29.0, -5.0, 0.0] dt = np.linspace(0.0, 86400 * 700, 5000) # time is seconds 

Answer 1

As you have shown, you can write this as a system of six first-order odes:

 x' = x2 y' = y2 z' = z2 x2' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * x y2' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * y z2' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * z 

You can save this as a vector:

 u = (x, y, z, x2, y2, z2) 

and thereby create a function that returns its derivative:

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

Given the initial state u0 = (x0, y0, z0, x20, y20, z20) and the variable for time t , this can be passed to scipy.integrate.odeint as such:

 u = odeint(deriv, u0, t) 

where u will be a list as above. Or you can unzip u from the beginning and ignore the values ​​for x2 , y2 and z2 (you must first transfer the output with .T )

 x, y, z, _, _, _ = odeint(deriv, u0, t).T