Deriving the Rocket Equation from First Principles

Wyatt Sheffield

The Tsiolkovsky Rocket Equation (hereafter referred to simply as “the rocket equation”) models the ultimate velocity of a rocket given its mass, the mass of fuel, and the velocity of the exhaust. While it may seem at first to be a straightforward relationship, the fundamental problem of rocketry is that a rocket must lift its own fuel. Say you want your 1000 kg rocket to accelerate to 1000 km/h, and you find that 200 kg of fuel has enough oomph to get 1000 kg up to that speed. So you fuel up your rocket, which now weighs 1200 kg, meaning it’s too heavy to get up to speed. So you have to add more fuel to lift the weight of the fuel. And again, and again and again…

The rocket equation prevents us from going through this iterative approach by giving an analytic solution to the problem. Risking accusations of spoilers, the rocket equation may be given as

Δ v = v e ln m 0 m f {\huge {\Delta v = v_e\ln\frac{m_0}{m_f}}}

in which

Where does this equation come from? I was familiar with it as a tool, but I had always had it handed down from on high without explanation. So, I decided to give it a shot myself, under the following constraints and assumptions:

Initial Analysis

We begin by considering a “rocket” that propels itself by launching a single large mass in the opposite direction of its travel. This rocket is essentially a cannon, using a cannonball as propellant. We can then consider a rocket that fires two successive cannonballs, each with half the mass of the cannonball in the first design. Then we may consider a 3-shot design, a 4-shot, etc.

Since the rocket starts from rest, and we are only allowed the law of conservation of momentum, we may write

∑ n p ⃗ n , i = ∑ m p ⃗ m , f = 0 \sum_n \vec{p}_{n,i} = \sum_m \vec{p}_{m,f} = 0

That is, the vector sum of all initial momenta at time t 0 t_0 must be zero and identical to the sum of all final momenta at time t 1 t_1 .

Case 1: A Single Propellant Ejection

The ejected propellant will have momentum m e v e m_ev_e that must be equal in magnitude to the momentum of the rocket m f Δ v m_f\Delta v . As we wish to solve for the final velocity of the rocket, we find

Δ v = m e v e m f . \Delta v = \frac{m_e v_e}{m_f}\text{.}

image/svg+xml

Case 2: Two sequential ejections

The first firing is calculated just as in the previous case, though we mustn’t forget that rocket still has extra mass from the remaining fuel.

1 2 m e v e = ( m f + 1 2 m e ) Δ v 1 Δ v 1 = m e v e 2 m f + m e \begin{align} \frac 1 2 m_e v_e &= (m_f + \frac 1 2 m_e) \Delta v_1 \\ \Delta v_1 &= \frac{m_e v_e}{2m_f + m_e} \end{align}

The second and final ejection occurs while the rocket is already in motion. We simply add Δ v 1 \Delta v_1 to the result of the second ejection.

1 2 m e v e = m f Δ v 2 Δ v 2 = m e v e 2 m f \begin{align} \frac 1 2 m_e v_e &= m_f \Delta v_2 \\ \Delta v_2 &= \frac{m_e v_e}{2m_f} \end{align}

Therefore

Δ v = Δ v 1 + Δ v 2 = m e v e 2 m f + m e + m e v e 2 m f \Delta v = \Delta v_1 +\Delta v_2 =\frac{m_e v_e}{2m_f + m_e} + \frac{m_e v_e}{2m_f}

Case 3

Following the same procedure, we can see that Δ v 1 = m e v e 3 m f + 2 m e \Delta v_1 = \frac{m_e v_e}{3m_f + 2m_e} and we can quickly conclude that the pattern extends to

Δ v = m e v e ( 1 3 m f + 2 m e + 1 3 m f + m e + 1 3 m f ) \Delta v = m_e v_e \left(\frac{1}{3m_f + 2m_e} + \frac{1}{3m_f + m_e} + \frac{1}{3m_f}\right)

Theorem

A rocket of mass m f m_f propelled by ejecting n n discreet “bullets” of mass 1 n m e \frac1n m_e at velocity v e v_e will attain a total change in velocity Δ v \Delta v given by Δ v = m e v e ∑ k = 1 n 1 n m f + ( n − k ) m e \big\displaystyle \Delta v = m_e v_e \sum_{k=1}^n \frac{1}{nm_f + (n-k)m_e}

Having derived this series formula, we now wish to find a closed-form analytic solution.

From Discrete to Continuous

We can imagine as we increase the number of fuel ejections that the process goes from a series of distinct firings to a continuous stream of mass being ejected from the rocket. Indeed this is physically true: as the number of projectiles increases the share of mass for each decreases until the projectiles are individual atoms. Thus, we must take the limit as n → ∞ n \to \infty to correctly model an actual rocket.

Δ v = lim n → ∞ m e v e ∑ k = 1 n 1 n m f + ( n − k ) m e \Delta v =\lim_{n\to\infty} m_e v_e \sum_{k=1}^n \frac{1}{nm_f + (n-k)m_e}

We can modify the equation slightly by letting x = m e m 0 x = \frac{m_e}{m_0} represent the proportion of total mass taken up by the fuel (note that this implies 0 ≤ x ≤ 1 0\le x \le 1 ). We can then use this substitution to rewrite the denominator

Δ v = lim n → ∞ m e v e ∑ k = 1 n 1 n m f + ( n − k ) m e = lim n → ∞ m e v e ∑ k = 1 n 1 n ( m 0 − m e ) + ( n − k ) m e = lim n → ∞ m e v e ∑ k = 1 n 1 m 0 [ n ( 1 − m e m 0 ) + ( n − k ) m e m 0 ] = lim n → ∞ m e v e ∑ k = 1 n 1 m 0 [ n ( 1 − x ) + ( n − k ) x ] = lim n → ∞ m e v e ∑ k = 1 n 1 m 0 [ n − n x + n x − k x ] = v e lim n → ∞ m e m 0 ∑ k = 1 n 1 n − k x = v e lim n → ∞ x ∑ k = 1 n 1 n − k x = v e lim n → ∞ ∑ k = 1 n x n − k x \begin{align} \Delta v&= \lim_{n\to\infty} m_e v_e \sum_{k=1}^n \frac{1}{nm_f + (n-k)m_e} \\ &= \lim_{n\to\infty} m_e v_e \sum_{k=1}^n \frac{1}{n(m_0 - m_e) + (n-k)m_e} \\ &= \lim_{n\to\infty} m_e v_e \sum_{k=1}^n \frac{1}{m_0\left[n\left(1-\frac{m_e}{m_0}\right) + (n-k)\frac{m_e}{m_0}\right]} \\ &= \lim_{n\to\infty} m_e v_e \sum_{k=1}^n \frac{1}{m_0\left[n(1-x) + (n-k)x\right]} \\ &= \lim_{n\to\infty} m_e v_e \sum_{k=1}^n \frac{1}{m_0\left[n-nx + nx - kx\right]} \\ &= v_e \lim_{n\to\infty} \frac{m_e}{m_0}\sum_{k=1}^n \frac{1}{n - kx} \\ &= v_e \lim_{n\to\infty} x \sum_{k=1}^n \frac{1}{n - kx} \\ &= v_e \lim_{n\to\infty} \sum_{k=1}^n \frac{x}{n - kx} \\ \end{align}

Though this expression is greatly simplified, we must contend with the unusual fact that the upper limit n n of the sum appears in the summation itself. We expect the sum to converge for some values of 0 ≤ x ≤ ∞ 0\le x \le \infty since for x = 0 x=0 , it is clear that lim n → ∞ ∑ k = 1 n x n − k x = 0 \displaystyle\lim_{n\to\infty} \sum_{k=1}^n \frac{x}{n - kx} = 0 and for x = 1 x=1 , the sum becomes identical to the harmonic series and diverges. Thus it is likely that there is some positive finite value that is attained on the interval 0 < x < 1 0\lt x \lt 1 . This by no means proves convergence, but merely gives evidence for the possibility as well as motivation to continue our analysis.

To evaluate such a “nearly-harmonic” series we may use the digamma function, ψ \psi , defined as the logarithmic derivative of the gamma function

ψ ( z ) = Γ ′ ( z ) Γ ( z ) = d d z ln Γ ( z ) \psi(z) = \frac{\Gamma'(z)}{\Gamma(z)} = \dd z \ln\Gamma(z)

The intuition for the usefulness of the digamma function is that the gamma function generalizes the factorial, and the log of a product (like a factorial) becomes a sum.

ln Γ ( x ) ∼ ln x ! = ln ( ∏ k = 1 x k ) = ∑ k = 1 x ln k = ∑ k = 1 x ln ( x − k − 1 ) \ln\Gamma(x) \sim \ln x! = \ln\left(\prod_{k=1}^x k\right)=\sum_{k=1}^x \ln k =\sum_{k=1}^x \ln(x-k-1)

Then by differentiating the above result, we have

ψ ( x ) = d d x ln Γ ( x ) ∼ d d x ∑ k = 1 x ln ( x − k − 1 ) = ∑ k = 1 x 1 x − k − 1 = ∑ k = 1 x 1 k = H x \psi(x)=\frac{\d}{\d x}\ln\Gamma(x) \sim \frac{\d}{\d x}\sum_{k=1}^x \ln(x-k-1) = \sum_{k=1}^x \frac{1}{x-k-1} = \sum_{k=1}^x \frac{1}{k} = H_x

Where H x H_x is the x th x\text{th} harmonic number. The actual value of ψ ( n ) \psi(n) for some natural number n n is ψ ( n ) = H n − 1 − γ \psi(n) = H_{n-1} - \gamma , supporting our intuition that the digamma function will be useful for computing a “nearly-harmonic” sum.

Theorem

Generalized recurrence relation for the digamma function

ψ ( z + N ) − ψ ( z ) = ∑ k = 0 N − 1 1 z + k \big\displaystyle {\psi(z+N) - \psi(z) = \sum_{k=0}^{N-1}\frac{1}{z + k}}

Proof: In order to use the digamma function to evaluate our sum, consider the defining recurrence relation of the gamma function Γ ( z + 1 ) = z Γ ( z ) \Gamma(z+1)=z\Gamma(z) . Taking the logarithmic derivative, we arrive at a recurrence relation for the digamma function

d d z ln Γ ( z + 1 ) = d d z ln ( z Γ ( z ) ) ψ ( z + 1 ) = d d z [ ln z + ln Γ ( z ) ] = 1 z + ψ ( z ) \begin{align} \dd{z} \ln\Gamma(z+1) &= \dd{z} \ln(z\Gamma(z)) \\ \psi(z+1) &= \dd{z} [\ln z + \ln\Gamma(z)] \\ &= \frac{1}{z} + \psi(z) \\ \end{align}

Now suppose that ψ ( z + N ) − ψ ( z ) = ∑ k = 0 N − 1 1 z + k \psi(z + N) - \psi(z) = \sum_{k=0}^{N-1}\frac{1}{z + k} . Then for N + 1 N+1 we have

ψ ( z + N + 1 ) = 1 z + N + ψ ( z + N ) = 1 z + N + ∑ k = 0 N − 1 1 z + k + ψ ( z ) = ∑ k = 0 N 1 z + k + ψ ( z ) □ \begin{align} \psi(z + N + 1) &= \frac{1}{z+N} + \psi(z+N) \\ &= \frac{1}{z+N} + \sum_{k=0}^{N-1}\frac{1}{z + k} + \psi(z) \\ &= \sum_{k=0}^{N}\frac{1}{z + k} + \psi(z) \quad\square \end{align}

We may now use this result to evaluate our rocket summation ∑ k = 1 n x n − k x \sum_{k=1}^n \frac{x}{n - kx} .

Δ v v e = ∑ k = 1 n x n − k x = − ∑ k = 1 n 1 − n x + k = − ∑ k = 0 n − 1 1 ( 1 − n x ) + k = ψ ( 1 − n x ) − ψ ( 1 + n − n x ) \frac{\Delta v}{v_e} = \sum_{k=1}^n \frac{x}{n - kx} = -\sum_{k=1}^n \frac{1}{-\frac{n}{x} + k} = -\sum_{k=0}^{n-1} \frac{1}{\left(1-\frac{n}{x}\right) + k} = \psi\left(1-\frac{n}{x}\right) - \psi\left(1+n-\frac{n}{x}\right)

Now we are ready to take the limit as n → ∞ n \to \infty using the fact that ψ ( z ) ∼ ln z − 1 2 z \psi(z)\sim\ln z - \frac{1}{2z} as z → ∞ z\to\infty .

Δ v v e = lim n → ∞ ∑ k = 1 n x n − k x = lim n → ∞ ψ ( 1 − n x ) − ψ ( 1 + n − n x ) = lim n → ∞ ln ( 1 − n x ) − 1 ( 1 − n x ) − [ ln ( 1 + n − n x ) − 1 ( 1 + n − n x ) ] = lim n → ∞ ln ( 1 − n x ) − ln ( 1 + n − n x ) = lim n → ∞ ln ( 1 − n x 1 + n − n x ) = lim n → ∞ ln ( 1 n − 1 x 1 n + 1 − 1 x ) = ln ( − 1 x 1 − 1 x ) = ln ( − 1 x − 1 ) = − ln ( 1 − x ) \begin{align} \frac{\Delta v}{v_e}&= \lim_{n\to\infty} \sum_{k=1}^n \frac{x}{n - kx} \\ &= \lim_{n\to\infty} \psi\left(1-\frac{n}{x}\right) - \psi\left(1+n-\frac{n}{x}\right) \\ &= \lim_{n\to\infty} \ln\left(1-\frac{n}{x}\right)-\cancel{\frac{1}{\left(1-\frac{n}{x}\right)}}-\left[ \ln\left(1+n-\frac{n}{x}\right) - \cancel{\frac{1}{\left(1+n-\frac{n}{x}\right)}}\right] \\ &= \lim_{n\to\infty} \ln\left(1-\frac{n}{x}\right)-\ln\left(1+n-\frac{n}{x}\right) \\ &= \lim_{n\to\infty} \ln\left(\frac{1-\frac{n}{x}}{1+n-\frac{n}{x}}\right) \\ &= \lim_{n\to\infty} \ln\left(\frac{\cancel{\frac1n}-\frac1x}{\cancel{\frac1n}+1-\frac1x}\right) = \ln\left(\frac{-\frac1x}{1-\frac1x}\right) = \ln\left(\frac{-1}{x-1}\right) \\ &= -\ln(1-x) \end{align}

Thus the rocket equation is Δ v = − v e ln ( 1 − m e m 0 ) \Delta v = -v_e\ln\left(1-\frac{m_e}{m_0}\right) . After checking my work against the accepted form of the equation (which can be seen in the introduction), I saw that mass ratio used by most authors was m 0 m f \frac{m_0}{m_f} . We can easily reframe our equation to match this form:

Δ v = − v e ln ( 1 − m e m 0 ) = − v e ln ( m 0 − m e m 0 ) = − v e ln ( m f m 0 ) \begin{align} \Delta v&= -v_e\ln\left(1-\frac{m_e}{m_0}\right) \\ &= -v_e\ln\left(\frac{m_0 - m_e}{m_0}\right) \\ &= -v_e\ln\left(\frac{m_f}{m_0}\right) \\ \end{align}

And we apply our last beautiful logarithm property to arrive precisely at our goal.

Δ v = v e ln ( m 0 m f ) \Delta v = v_e\ln\left(\frac{m_0}{m_f}\right)

Closing thoughts

Other derivations of the rocket equation are certainly more elegant as they can take advantage of other physical laws. However, this nearly brute-force approach turned out to be deeply enlightening - I never would have expected the gamma function to come into play, for instance. Furthermore, I had never before had a chance to explore the digamma function. Without solving this problem, I doubt I would have had a chance to learn about it.

Solve hard problems in the wrong way. You’ll always learn something.

Scans of Relevant pages from my Notebook

Page 1 of my original notes on the topic Page 2 of my original notes on the topic Page 3 of my original notes on the topic Page 4 of my original notes on the topic
Tags
math research physics calculus rockets infinite series special functions