# Robust FindRoot search

I have the following simulation that models a rocket landing (flat Earth, no atmosphere, constant gravitational acceleration), which uses root finding to figure out when to start and stop the engines from firing.

Remove["Global*"]
g = 9.81; (*Gravitational acceleration*)
mprop = (409500)*0.15 ;(*15% of propellant mass left after MECO*)
mstruct = 22200; (*Structural mass*)
mtotal = mprop + mstruct;(*Total mass*)
T[n_] := 845000 n;(*Thrust from n engines*)
Isp = 282; (*Specific impulse*)
tmax = 1000;(*Maximum value for t*)

theta = Pi;(*Initial thrust angle in radians*)

Solution = ParametricNDSolve[{
x''[t] == If[t < tburn1end, T/m[t] Cos[theta], If[tburn2start < t, -(T/m[t]) Cos[gamma[t]], If[m[t] < mstruct, 0, 0]]],
y''[t] == If[t < tburn1end, T/m[t] Sin[theta], If[tburn2start < t, -(T/m[t]) Sin[gamma[t]], If[m[t] < mstruct, 0, 0]]] - g,
m'[t] == If[t < tburn1end, -(T/(g Isp)), If[tburn2start < t, -(T/(g Isp)), If[m[t] < mstruct, 0, 0]]], gamma'[t] == D[ArcTan[x'[t], y'[t]], t],
x == 30000, y == 90000, x' == 1900 Cos[70 Degree], y' == 1900 Sin[70 Degree], m == mtotal, gamma == 70 Degree},
{x[t], y[t], x'[t], y'[t], x''[t], y''[t], m[t], gamma[t]}, {t, 0, tmax}, {tburn1end, tburn2start},
MaxSteps -> 1000000(*, Method\[Rule]"StiffnessSwitching"*)]

rootxy = FindRoot[{x[t][tburn1end, tburn2start] == 0, y[t][tburn1end, tburn2start] == 0, y'[t][tburn1end, tburn2start] == 0} /.Solution, {{t, 400}, {tburn1end, 25}, {tburn2start, 300}}]
tburn1end = rootxy[][]
tburn2start = rootxy[][]
tburn2end = rootxy[][]

ParametricPlot[Evaluate[{x[t][tburn1end, tburn2start], y[t][tburn1end, tburn2start]} /. Solution], {t, 0, tburn2end}, AxesLabel -> {x, y}, PlotRange -> {{-20000, 200000}, {-10000, 300000}}, PlotStyle -> Automatic, ImageSize -> Large]
Animate[ParametricPlot[Evaluate[{x[t][tburn1end, tburn2start], y[t][tburn1end, tburn2start]} /. Solution], {t, 0, a}, AxesLabel -> {x, y},PlotRange -> {{-20000, 200000}, {-10000, 300000}}, PlotStyle -> Automatic, ImageSize -> Large], {a, 0, tburn2end}, AnimationRate -> 5, AnimationRepetitions -> 1]

Plot[{x[t][tburn1end, tburn2start]} /. Solution, {t, 0, tburn2end},AxesLabel -> {t, x}]
Plot[{y[t][tburn1end, tburn2start]} /. Solution, {t, 0, tburn2end},AxesLabel -> {t, y}]
Plot[{x'[t][tburn1end, tburn2start]} /. Solution, {t, 0, tburn2end}, AxesLabel -> {t, vx}]
Plot[{y'[t][tburn1end, tburn2start]} /. Solution, {t, 0, tburn2end}, AxesLabel -> {t, vy}]


However, I've noticed that FindRoot isn't very robust, and the guessed values have to be fairly close to the actual roots for a solution to be found. It would be great if I was able to set my search limits from {t,0,tmax} and use something like Brent's method which, as far as I know, is guaranteed to find a solution if it exists within the search bounds due to its use of the bisection method. But unfortunately Mathmatica says it only works with univariate functions. As such, I'm wondering if there is another way to be able to find the required engine start and stop times that is more reliable than what I'm currently doing. I've had a go at changing various settings in ParametricNDSolve and FindRoot, such as PrecisionGoal and AccuracyGoal, but nothing has seemed to work so far. Any help would be great, thanks!

EDIT: In response to Michael E2, I tried using NMinimize[] instead of FindRoot[], but my results were even worse. Here's the code I used in case I did something wrong which may have produced the bad results:

QuadraticMin[t0_?NumericQ, tburn1end_?NumericQ, tburn2start_?NumericQ] := Module[{tem = (x[t][tburn1end, tburn2start])^2 + (y[t][tburn1end,tburn2start])^2 + (y'[t][tburn1end, tburn2start])^2 /.Solution}, tem /. t -> t0];

Constraints[t0_?NumericQ, tburn1end_?NumericQ, tburn2start_?NumericQ] := Module[{tem = (y[t][tburn1end, tburn2start]) /. Solution}, tem /. t -> t0];

minimize = NMinimize[{QuadraticMin[t, tburn1end, tburn2start],Constraints[t, tburn1end, tburn2start] > 0, 0 < tburn1end < tmax, 0 < tburn2start < tmax, 0 < t < tmax}, {{t, 0, tmax}, {tburn1end, 0, tmax}, {tburn2start, 0, tmax}}]


The results were as follows:

{9.00056*10^9, {t -> 7.82522*10^-9, tburn1end -> 406.714, tburn2start -> 1.6873*10^-8}}

• You might want to use WhenEvent[] instead to determine when your components become zero. – J. M. is away Jul 27 '16 at 10:09
• Also, take a look at mathematica.stackexchange.com/questions/91784/… – yohbs Jul 27 '16 at 10:50
• @J.M. I'm a little confused as to how I'd use WhenEvent[] to find the values for the parameters tburn1end and tburn2start that will give me values of x[t] == 0 && y[t] == 0 && y'[t] == 0. I thought WhenEvent[] was used to set values when certain events occur and could not be used to find parameters that produce desired results. – InquisitiveInquirer Jul 27 '16 at 12:57
• Have you tried NMinimize[] on $x^2+y^2+(y')^2$ (or perhaps FindMinimum[])? – Michael E2 Jul 27 '16 at 12:59
• @MichaelE2, unfortunately NMinimize[] seems to fare even worse. What is the difference between NMinimize[] and FindMinimum[]`? – InquisitiveInquirer Jul 27 '16 at 13:33