# Problem using ArcTan inside NDSolve

I'm messing around a bit with using "thrust segments" on a point-mass in NDSolve and have run into a strange problem to do with using ArcTan inside NDSolve. The code below applies thrust to a point-mass for 30 seconds, and then lets it move in free-fall from then onward.

Remove["Global*"]
g = 9.81; (*Gravitational acceleration*)
m0 = 50000; (*Initial mass*)
T = 1200000;(*Thrust*)
Isp = 300; (*Specific impulse*)
tmax = 1000;(*Maximum value for t*)
theta = 177;(*Thrust angle*)

Solution = NDSolve[{
x''[t] == If[t < 30, T/m[t] Cos[theta Degree], 0],
y''[t] == If[t < 30, T/m[t] Sin[theta Degree], 0] - g,
m'[t] == If[t < 30, -(T/(g Isp)), 0],
gamma[t] == ArcTan[x'[t], y'[t]],
x[0] == 100000, y[0] == 100000, x'[0] == 1900 Cos[70 Degree],
y'[0] == 1900 Sin[70 Degree], m[0] == m0}, {x[t], y[t], x'[t],
y'[t], m[t], gamma[t]}, {t, 0, tmax}, MaxSteps -> 1000000]

gammatable = Table[ArcTan[x'[t], y'[t]]*180/\[Pi] /. Solution, {t, 0, tmax}]

ParametricPlot[Evaluate[{x[t], y[t]} /. Solution], {t, 0, tmax},
AxesLabel -> {x, y}, PlotRange -> {{0, 900000}, {0, 300000}},
PlotStyle -> Automatic, ImageSize -> Large]
Animate[ParametricPlot[Evaluate[{x[t], y[t]} /. Solution], {t, 0, a},
AxesLabel -> {x, y}, PlotRange -> {{0, 900000}, {0, 300000}},
PlotStyle -> Automatic, ImageSize -> Large], {a, 0, tmax},
AnimationRate -> 5, AnimationRepetitions -> 1]


As can be seen, a term called gamma is calculated inside NDSolve, which is just the flight-path angle of the point mass, gamma[t] == ArcTan[x'[t], y'[t]]. It seems that when the y-velocity of the particle becomes zero (at the top of the free-fall arc), the calculation of gamma goes haywire and the subsequent plot and animation "break". Looking at the values that gamma attains at the point where the y-velocity is zero (produced using gammatable), we see that the value of ArcTan is very close to zero. However, from then onward the value stays close to zero, even though the particle gains negative y-velocity. As such, I'm wondering if there's a way to prevent this from happening so that an accurate value for gamma is given at each point in time from t = 0 to t = tmax. I tried using the following gamma[t] == If[y'[t] < Abs[0.01], \[Pi], ArcTan[x'[t], y'[t]]], which seemed to fix the problem. However this quick-fix is less than ideal and I was hoping there there is a cleaner way using no tricks to fix the problem. Any help would be great, thanks!

• Are you sure you want to be using "remove"? The basic methods for stiffness switching unfortunately don't work here. Jul 24 '16 at 20:18

One way to avoid the discontinuity of the branch cut in ArcTan is to integrate its derivative. This also converts the DAE, which is restricted to use the IDA method, to an ODE. This allows greater flexibility in integrating the system. In this case, the default settings seem to work well as is.

g = 9.81;      (*Gravitational acceleration*)
m0 = 50000;    (*Initial mass*)
T = 1200000;   (*Thrust*)
Isp = 300;     (*Specific impulse*)
tmax = 1000;   (*Maximum value for t*)
theta = 177;   (*Thrust angle*)
{Solution} =
NDSolve[{     (* Piecewise[] is a better choice than If[] for functions/equations *)
x''[t] == Piecewise[{{T/m[t] Cos[theta Degree], t < 30}}],
y''[t] == Piecewise[{{T/m[t] Sin[theta Degree], t < 30}}] - g,
m'[t] == Piecewise[{{-(T/(g Isp)), t < 30}}],
gamma'[t] == D[ArcTan[x'[t], y'[t]], t],        (* <-- N.B. *)
x[0] == 100000, y[0] == 100000,
x'[0] == 1900 Cos[70 Degree], y'[0] == 1900 Sin[70 Degree],
m[0] == m0,
gamma[0] == 70 Degree},                         (* <-- N.B. *)
{x, y, m, gamma}, {t, 0, tmax}]

gammatable = Table[ArcTan[x'[t], y'[t]] /. Solution, {t, 0, tmax}];


Plot of gamma[t]:

Plot[Evaluate[gamma[t] /. Solution], {t, 0, tmax}]


Compare with the original idea for gamma[t] as ArcTan[x'[t], y'[t]]:

ListPlot@gammatable


Aside from a glitch at the discontinuity in the ODE at t == 30, there is a discontinuity in ArcTan[x'[t], y'[t]] around t == 186, which is where the OP's original NDSolve[] ran into a stiffness/singularity problem (NDSolve::ndsz).

ParametricPlot[Evaluate[{x[t], y[t]} /. Solution], {t, 0, tmax},
AxesLabel -> {x, y}, PlotRange -> {{0, 300000}, {0, 300000}}]


• Nice ... Ultimately it looks like the same answer. Jul 25 '16 at 5:08
• @Young They shouldn't be too far apart. It's hard to say what is the net effect of reducing PrecisionGoal and raising AccuracyGoal (in your method). Theoretically, PrecisionGoal -> 1 would allow as much as a ~10% error in a local step, but the two solutions do not appear that far off from each other. Jul 25 '16 at 5:19
• Very nice little "trick", works like a charm. Jul 25 '16 at 10:12
• As an aside, do you know why Piecewise is better than using If in this situation? Jul 25 '16 at 11:53
• @user7388 Piecewise is intended for constructing functions and If is for programming (control flow). In particular, solvers treat the boundaries of cases in Piecewise as potential discontinuities and adapt the solving method to the discontinuities. This is important because error estimators generally assume the function has continuous derivatives up to some order (depending on the order of the method). In recent versions, some/all solvers seem to handle some If statements, so my advice is perhaps less important now. It's more of a precaution. Jul 25 '16 at 14:40

One of the methods I've used with stiff systems is to reduce the PrecisionGoal (and sometimes AccuracyGoal). Check the result for reasonableness.

Solution = NDSolve[{
x''[t] == If[t < 30, T/m[t] Cos[theta Degree], 0],
y''[t] == If[t < 30, T/m[t] Sin[theta Degree], 0] - g,
m'[t] == If[t < 30, -(T/(g Isp)), 0],
gamma[t] == ArcTan[x'[t], y'[t]],
x[0] == 100000, y[0] == 100000, x'[0] == 1900 Cos[70 Degree],
y'[0] == 1900 Sin[70 Degree], m[0] == m0},
{x[t], y[t], x'[t], y'[t], m[t], gamma[t]}, {t, 0, tmax},
MaxSteps -> Infinity, PrecisionGoal -> 1, AccuracyGoal -> 10];


arc = Table[
ParametricPlot[Evaluate[{x[t], y[t]} /. Solution], {t, 0, a},
AxesLabel -> {x, y}, PlotRange -> {{0, 300000}, {0, 300000}},
PlotLabel -> a, PlotStyle -> Automatic, ImageSize -> Large],
{a, 1, 400}];
Export["arc.gif", arc]


Note:
The above result agrees with gamma[t] == If[y'[t] < Abs[0.01], \[Pi], ArcTan[x'[t], y'[t]]]`