8
$\begingroup$

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
(*Radius of Earth*)
r = AstronomicalData["Earth", "Radius"]
(*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.

Improve this question
$\endgroup$
3
  • $\begingroup$ If you supply the conditional statements you want added, someone will show you how. $\endgroup$
    – Nasser
    Commented Aug 22, 2013 at 21:51
  • $\begingroup$ I'm looking to do something similar to the following: <code> if(t==1000){ T==0 } </code> $\endgroup$
    – InquisitiveInquirer
    Commented Aug 22, 2013 at 21:56
  • 2
    $\begingroup$ 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. $\endgroup$
    – Nasser
    Commented Aug 22, 2013 at 22:01

1 Answer 1

Reset to default
8
$\begingroup$

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]}]

enter image description here

Improve this answer
$\endgroup$
8
  • $\begingroup$ Thanks Kuba, this looks perfect! $\endgroup$
    – InquisitiveInquirer
    Commented Aug 23, 2013 at 6:47
  • $\begingroup$ @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. $\endgroup$
    – Kuba
    Commented Aug 23, 2013 at 6:53
  • $\begingroup$ 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? $\endgroup$
    – InquisitiveInquirer
    Commented Aug 23, 2013 at 7:13
  • $\begingroup$ @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. $\endgroup$
    – Kuba
    Commented Aug 23, 2013 at 7:25
  • $\begingroup$ 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]? $\endgroup$
    – InquisitiveInquirer
    Commented Aug 23, 2013 at 11:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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