Solving a second order non-linear ODE (rocket height function) $h''(t)=\frac{F}{m+k-r\cdot t}-\frac{G\cdot M}{(R+h)^{2}}$

Question

I've talked to a few people, and most of them think that it is unlikely that there is an exact solution to this differential equation.

Here is the equation:

$$h''(t)=\frac{F}{m+k-r\cdot t}-\frac{G\cdot M}{(R+h)^{2}}$$

$F,m,k,r,G,M,R$ are all constants. I'm trying to find solve for $h(t)$, which represents the altitude of a rocket. If it makes it any simpler, I'm only trying to find $h$ when $t=\frac{k}{r}$. I'm also wondering whether an approximation can be obtained, if there are no solutions.

Background (Skip if you want)

Imagine a single stage rocket with a dry mass of $m$ which produces a fixed thrust of $F$ Newtons. It carries with it $k$ kilograms of fuel and uses $r$ kilograms of fuel per second. Using Newton's 2nd Law $F=m\cdot a$, and assuming the rocket simply flies straight up, the acceleration of the rocket is:

$$a(t)=\frac{F}{m+k-r\cdot t}-g$$

where $t$ is the elapsed time in seconds and $g$ is the downwards acceleration due to gravity.

I would like to find the height of the rocket when its fuel is expended (when $t=\frac{k}{r}$). However, I also want the equation to account for the decreasing gravitational force as the rocket increases in altitude:

$$a(t)=\frac{F}{m+k-r\cdot t}-\frac{G\cdot M}{(R+h)^{2}}$$

$G$ is the gravitational constant, $M$ and $R$ is the mass and radius of the earth respectively, and $h$ is the altitude of the rocket. Since acceleration is defined as the second derivative of displacement, or in this case altitude or height, the resulting equation is:

$$h''(t)=\frac{F}{m+k-r\cdot t}-\frac{G\cdot M}{(R+h)^{2}}$$

Answer 1

At first we solve by numeric method.(Analytically can't be solved. Maple and Mathematica can't solve)

Mathematica code:

G = 6.67408*10^-11;
M = 5.9722*10^24;
R = 6378.14*1000;
F = 68000;
m = 2340;
k = 3000;
r = 68;
t0 = k/r // N

when $t=k/r$ gives: 44.11764706 second. $$h''(t)=\frac{F}{k+m-r t}-\frac{G M}{(h(t)+R)^2}$$

sol = NDSolve[{h''[t] == F/(m + k - r *t) - (G*M)/(R + h[t])^2, 
h[0] == 0, h'[0] == 0}, h, {t, 0, t0}];
((h[t] /. sol) /. t -> t0)

putting numeric values to solution:

Rocket will lift at the height: 6192.405405 meters (6.2 km).

---

At second we will a approximation at point h=0 expanding with series:

$$\frac{G M}{(h(t)+R)^2}=\frac{G M}{R^2}-\frac{2 (G M) h(t)}{R^3}+O\left(h(t)^2\right)$$

putting to equation: $$h''(t)=\frac{F}{k+m-r t}+\frac{2 G M h(t)}{R^3}-\frac{G M}{R^2}$$

With help CAS we can solve:

$h(t)=-\frac{R e^{-\frac{2 \sqrt{2} \sqrt{G} \sqrt{M} (k+m+r t)}{r R^{3/2}}} \left(\sqrt{2} F \sqrt{R} \text{Ei}\left(-\frac{\sqrt{2} \sqrt{G} (k+m) \sqrt{M}}{r R^{3/2}}\right) e^{\frac{\sqrt{2} \sqrt{G} \sqrt{M} (3 k+3 m+r t)}{r R^{3/2}}}-\sqrt{2} F \sqrt{R} \text{Ei}\left(\frac{\sqrt{2} \sqrt{G} (k+m) \sqrt{M}}{r R^{3/2}}\right) e^{\frac{\sqrt{2} \sqrt{G} \sqrt{M} (k+m+3 r t)}{r R^{3/2}}}-\sqrt{2} F \sqrt{R} e^{\frac{\sqrt{2} \sqrt{G} \sqrt{M} (3 k+3 m+r t)}{r R^{3/2}}} \text{Ei}\left(-\frac{\sqrt{2} \sqrt{G} \sqrt{M} (k+m-r t)}{r R^{3/2}}\right)+\sqrt{2} F \sqrt{R} e^{\frac{\sqrt{2} \sqrt{G} \sqrt{M} (k+m+3 r t)}{r R^{3/2}}} \text{Ei}\left(\frac{\sqrt{2} \sqrt{G} \sqrt{M} (k+m-r t)}{r R^{3/2}}\right)-2 \sqrt{G} \sqrt{M} r e^{\frac{2 \sqrt{2} \sqrt{G} \sqrt{M} (k+m+r t)}{r R^{3/2}}}+\sqrt{G} \sqrt{M} r e^{\frac{\sqrt{2} \sqrt{G} \sqrt{M} (2 k+2 m+r t)}{r R^{3/2}}}+\sqrt{G} \sqrt{M} r e^{\frac{\sqrt{2} \sqrt{G} \sqrt{M} (2 k+2 m+3 r t)}{r R^{3/2}}}\right)}{4 \sqrt{G} \sqrt{M} r}$

 sol2 = DSolve[{h''[t] == 
 F/(m + k - r *t) - (G M)/R^2 + (2 (G M) h[t])/R^3, h[0] == 0, 
 h'[0] == 0}, h[t], t] // FullSimplify
 h[t] /. sol2 /. t -> t0

 (* {{h[t] -> -(1/(4 Sqrt[G] Sqrt[M] r))
E^(-((2 Sqrt[2] Sqrt[G] Sqrt[M] (k + m + r t))/(r R^(3/2))))
  R (E^((Sqrt[2] Sqrt[G] Sqrt[M] (2 (k + m) + r t))/(
    r R^(3/2))) (-1 + E^((Sqrt[2] Sqrt[G] Sqrt[M] t)/R^(
      3/2)))^2 Sqrt[G] Sqrt[M] r + 
   Sqrt[2] F Sqrt[
    R] (E^((Sqrt[2] Sqrt[G] Sqrt[M] (3 (k + m) + r t))/(
       r R^(3/2))) (ExpIntegralEi[-((
           Sqrt[2] Sqrt[G] (k + m) Sqrt[M])/(r R^(3/2)))] - 
         ExpIntegralEi[-((Sqrt[2] Sqrt[G] Sqrt[M] (k + m - r t))/(
           r R^(3/2)))]) + 
      E^((Sqrt[2] Sqrt[G] Sqrt[M] (k + m + 3 r t))/(
       r R^(3/2))) (-ExpIntegralEi[(
           Sqrt[2] Sqrt[G] (k + m) Sqrt[M])/(r R^(3/2))] + 
         ExpIntegralEi[(Sqrt[2] Sqrt[G] Sqrt[M] (k + m - r t))/(
          r R^(3/2))])))}}*)

putting numeric values to solution:

Rocket will lift at the height: 6192.406503 meters.

The difference form numerics and analytically is very small about 1 millimeter!

---

Another approximation.See-> Altitude:

$$h''(t)=\frac{F}{k+m-r t}-g$$ where $g$ is standard acceleration of free fall.

g = (G*M)/R^2
(* 9.798004977 *)

and solution:

$$h(t)=\frac{r t (2 F-g r t)+2 F (k+m-r t) (\log (k+m-r t)-\log (k+m))}{2 r^2}$$

sol3 = DSolve[{h''[t] == F/(m + k - r *t) - g, h[0] == 0, h'[0] == 0},
h[t], t] // FullSimplify
(h[t] /. sol3 /. t -> t0)

(* {{h[t] -> (r t (2 F - g r t) + 
2 F (k + m - r t) (-Log[k + m] + Log[k + m - r t]))/(2 r^2)}}*)

Rocket will lift at the height: 6190.114097 meters.

The difference form numerics and analytically is about 2 meters!