g = 0.0098; R = 6350; τ= 86; p = 0.75;a = 0.0446; s = 86.5 Degree; T = 3960;(* g = 0.0098; R = 6350; τ = 87; p = 0.75; a= 0.0446; s = 45 Degree; T = 5100; *)fr[t_, τ_] := a*Sin[s]*HeavisideTheta[τ - t]ft[t_, τ_] := a*Cos[s]*r[t]*HeavisideTheta[τ - t]fm[t_, τ_] := -p/(τ - p*t)*HeavisideTheta[τ - t]sol = NDSolve[{r''[t] == -fm[t, τ]*r'[t] + h[t]^2/r[t]^3 -g R^2/r[t]^2 + fr[t, τ], r[0] == R, r'[0] == 0,h'[t] == -fm[t, τ]*h[t] + ft[t, τ], h[0] == 0}, {r, h},{t, 0, T}]f[t_] := r[t] /. sol[[1, 1]]d[t_] := h[t] /. sol[[1, 2]]Plot[{6350, f[t]}, {t, 0, T}]Plot[{f[t + 1] - f[t], d[t]*R/f[t]^2, d[t]/f[t]},{t, 0, T}, PlotRange -> {-4, 8}]tyx[T_] := {T, f[T] - R, NIntegrate[R d[t]/f[t]^2 , {t, 0, T}]}v[T_] := Sqrt[(f[T] - f[T - 1])^2 + (R d[T]/f[T]^2)^2]{tyx[T], v[T]/0.340}g0 = ParametricPlot[{NIntegrate[R d[t]/f[t]^2 , {t, 0, tt}],f[tt] - R}, {tt, 0, T}]
燃焼時間τが86秒,燃料重量比が0.75,燃焼加速度が 0.0446km/s^2,投射角が86.5度である。これで飛行時間を与えると,飛行距離と最高高度と落下時速度が概ね再現できる。万有引力は距離の2乗に反比例し,コリオリ力や空気抵抗は無視,地球の形を考慮して2次元極座標で質量が変化する1段ロケットの微分方程式を解いている。
燃焼時間τを87秒にして,投射角を45度にすれば,飛行時間が5100秒で飛行距離は14,600km,落下時速度はマッハ24になる。このときペイロードはほとんど変化させていない。燃料重量比は同じで燃焼時間を1秒のばしただけだ。