The mass of a rocket as a function of time is: $$ m = m_0 - \dot m t$$ Where $m_0 $ is the initial mass of the rocket and $\dot m$ is the mass flow rate. At a certain time all the fuel has been used, so from that time onwards $m = m_f$, the mass of the rocket without fuel.

The speed is then defined by the following equation, where $V_e$ is the exhaust velocity and $g$ is constant gravity. $$v = V_e log({m_o \over m_o - \dot mt}) - gt$$

To obtain position let $\dot m$ be constant, the equation can be written as

$$v = V_e log({m_o \over m_o - kt}) - gt$$ Integration of that expression lets to the following, which corresponds with the expression contained in Orbital Mechanics by Curtis.

$$h = \frac{V_e}{\dot m}[(m_0 - \dot m t)log({m_o - \dot mt\over m_o })+ \dot mt] - \frac{1}{2}g t^2 $$


However, the plot of $h$ (Red) is different for what I get for evaluating $v$ and then doing a numerical integration (Blue). Both plots are similar until the time where the fuel runs out, then the difference is easily visible.

The analytic plot sounds more reasonable to me since a rocket needs multiple stages to achieve a high altitude, but the numerical method is harder to get wrong and is necessary to consider other factors such as drag. Can someone point me in the right direction? Am I doing something wrong?

Here is the code I used in MATLAB. The values used are an approximation to the Saturn V, just considering one stage.

t = 0:1:1600;

mass_inital = 2970.000;
mass_final = 731.000;
flow = 13.120;
exhaust = 2.500;
g = 9.81/1000;

mass = mass_inital - flow .* t;
mass(mass< mass_final) = mass_final;

v = exhaust .* log(mass_inital./mass) - (g .* t);

c = cumtrapz(v);
h = (exhaust ./ flow).* (mass.*log(mass./ mass_inital)+ flow .* t) - (0.5 .*g .*(t.*t));

hold on
ylim([0 max(c)]);

enter image description here


closed as off-topic by Kyle Kanos, user191954, John Rennie, Jon Custer, ZeroTheHero Oct 17 '18 at 0:15

  • This question does not appear to be about physics within the scope defined in the help center.
If this question can be reworded to fit the rules in the help center, please edit the question.

  • 4
    $\begingroup$ I'm voting to close this question as off-topic because it's about debugging code and not physics. $\endgroup$ – Kyle Kanos Oct 15 '18 at 21:36
  • $\begingroup$ "the numerical method is harder to get wrong" - that is one of the funnier things I've read in a while! People get numerical methods wrong all the time... $\endgroup$ – Jon Custer Oct 16 '18 at 19:19

Your fourth equation (for h) in the background section indeed seems to be the analytical integral of v. However, this is only true for the limited duration until fuel is exhausted.

From there on, the equation for v should be modified and this would modify the integral for h.

Most likely the last term in the square brackets for h ($\dot{m}t$) is incorrect when fuel is exhausted and is imposing extra height.

  • 1
    $\begingroup$ Yes! Adding the following equations after the time the fuel runs out (tbo) yielded the curve obtained in numerical integration. Thanks; v = (vbo - g.*(t-tbo)); h = (hbo + vbo.*(t-tbo)- (0.5 .* g .*power((t-tbo),2))); $\endgroup$ – John Oct 16 '18 at 0:58

Not the answer you're looking for? Browse other questions tagged or ask your own question.