# Modelling a Rocket Launch using NDSolve

I'm trying to model a rocket launch with Mathematica but I've run into a little problem since I don't know how to turn the thrust off. I'm using Newton's Law of Universal Gravitation plus an added term for the rocket's thrust and trying to find its position as a function of time. Here is the code so far (it is at a very basic stage at the moment with constant mass and no drag added, plus the thrust never stops, which is my biggest concern right now):

(*Gravitational Constant*)
G = 6.672*10^-11
(*Mass of Earth and rocket*)
M = AstronomicalData["Earth", "Mass"]
m = 2800000
(*Rocket thrust*)
T = 34020000
(*Numerical solution modelling the gravitation interation between the \
Earth and a launching rocket*)
(*NOTE: Rocket mass will change over time; also, add in drag*)
soln = NDSolve[{
x''[t] == -((G M x[t])/(x[t]^2 + y[t]^2 + z[t]^2)^(3/2)),
y''[t] == -((G M y[t])/(x[t]^2 + y[t]^2 + z[t]^2)^(3/2)) + 0.25 T/m,
z''[t] == -((G M z[t])/(x[t]^2 + y[t]^2 + z[t]^2)^(3/2)) + 0.75 T/m,
x[0] == 0, y[0] == 0, z[0] == r, x'[0] == 0, y'[0] == 0,
z'[0] == 0}, {x[t], y[t], z[t]}, {t, 0, 20000},
MaxSteps -> 10000000, Method -> "StiffnessSwitching"]

RocketPlot =
ParametricPlot3D[{x[t], y[t], z[t]} /. soln, {t, 0, 10000000}, AxesLabel -> {x, y, z}]


Does anyone know if conditional statements can be used inside NDSolve so that the thrust can be stopped at a certain position or time? Any help would be appreciated, cheers guys.

## migrated from mathematica.meta.stackexchange.comAug 22 '13 at 21:34

This question came from our discussion, support, and feature requests site for users of Wolfram Mathematica.

• If you supply the conditional statements you want added, someone will show you how. – Nasser Aug 22 '13 at 21:51
• I'm looking to do something similar to the following: <code> if(t==1000){ T==0 } </code> – InquisitiveInquirer Aug 22 '13 at 21:56
• So you want to integrate until t=1000, then want T=0, then? and after that, want to continue integrating? Then why not simply integrate up to t=1000 to start with, then use the state at the end, to start new NDSolve with now T=0 set into the equations, but using that state as initial conditions for next stage? It is also not always a good idea to integrate for very long span. I found it better to integrate over smaller time periods. This depends on your system ofcourse and how stiff it is or not. – Nasser Aug 22 '13 at 22:01

Yes you can, for example:

thrust[t_, t0_: 1000] := 34020.000 UnitStep[t0 - t]

end = 10000

soln = Table[
NDSolve[{
x''[t] == -((G M x[t])/Norm[{x[t], y[t], z[t]}]^3),
y''[t] == -((G M y[t])/Norm[{x[t], y[t], z[t]}]^3) + 0.25 thrust[t, t0]/m,
z''[t] == -((G M z[t])/Norm[{x[t], y[t], z[t]}]^3) + 0.75 thrust[t, t0]/m,
x[0] == 0, y[0] == 0, z[0] == r, x'[0] == 0, y'[0] == 0, z'[0] == 0},
{x[t], y[t], z[t]}, {t, 0, end}, MaxSteps -> 10000000,
Method -> "StiffnessSwitching"][[1]]
, {t0, 1000, 9000, 2000}]

ParametricPlot3D[Evaluate[{x[t], y[t], z[t]} /. soln], {t, 0, end},
AxesLabel -> {x, y, z}, AspectRatio -> 1, BoxRatios -> 1,
PlotStyle -> Automatic, BaseStyle -> {Thickness[.005]}]


• Thanks Kuba, this looks perfect! – InquisitiveInquirer Aug 23 '13 at 6:47
• @user7388 I hope it will help. Also, you can use If or something similar, it's just one example. But what Nasser has written is true, if thrust is going to be disabled most of the time it may be good idea to split integration to not check this condidtion when not needed. Good luck. – Kuba Aug 23 '13 at 6:53
• Thanks, I'll look into using an if statement and splitting integration up as well (unless someone would be willing to give me a hint as to what it would look like :P). Could you explain the thrust function a little further as I'm a bit confused by it. It looks like you're making the thrust function take two arguments, t and t0, and assigning 1000 to t0 in the definition. The part that most confuses me though is the use of the UnitStep function and [t0-t]; what's happening there? – InquisitiveInquirer Aug 23 '13 at 7:13
• @user7388 Take a look at UnitStep. t0 is just a shift so the jump is there and there is - t so it is 1 for t < t0 not for t > t0. At the end is T or 0. – Kuba Aug 23 '13 at 7:25
• I thought this might be a good time to ask a certain question that's been bothering me for a while: Why, when breaking up F=GMm/r^2 into its 3 vector equations, do we have to multiply each component equation by a unit vector? I know that unit vectors have a magnitude of 1, but for example if I take x''[t] == -GM/(Norm[x[t],y[t],z[t]]) surely this would be sufficient to find the acceleration in the x-direction since it is implied by x''[t]? – InquisitiveInquirer Aug 23 '13 at 11:28