Skip to content

Instantly share code, notes, and snippets.

@rikusalminen
Created July 31, 2014 09:28
Show Gist options
  • Save rikusalminen/6d6bb834d48b9664b38d to your computer and use it in GitHub Desktop.
Save rikusalminen/6d6bb834d48b9664b38d to your computer and use it in GitHub Desktop.
Solving the rocket equation in Python
def numerical():
from scipy.integrate import odeint
from numpy import array
# [position, mass, velocity, mass flow]
def acceleration(y, t, v_exhaust=1.0):
return array([y[2], y[3], v_exhaust * y[3]/y[1], 0.0])
initial_state = array([0.0, 2.0, 0.0, -1.0])
solution = odeint(acceleration, initial_state, [i/10.0 for i in range(11)], args=(1.0,))
print solution
def symbolic():
from sympy import symbols, integrate
m, n, t, v_exhaust = symbols('m n t v_exhaust')
acceleration = -v_exhaust * n / (m - n*t)
velocity = integrate(acceleration, (t, 0, t))
position = integrate(velocity, (t, 0, t))
params = {m: 2.0, n: 1.0, v_exhaust: 1.0}
v = velocity.subs(params)
p = position.subs(params)
solution = [(p.subs({t: i/10.0}).evalf(), v.subs({t: i/10.0}).evalf()) for i in range(11)]
print '\n'.join(map(str, solution))
def analytical():
from math import log
v_exhaust = 1.0
m = 2.0
n = 1.0
def velocity(t):
return v_exhaust * log((m - n*t) / m)
def position(t):
return -(v_exhaust / n) * ((m-n*t) * log((m-n*t) / m) + n*t)
solution = ((position(i/10.0), velocity(i/10.0)) for i in range(11))
print '\n'.join(map(str, solution))
numerical()
symbolic()
analytical()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment