# NDSolve how to monitor shooting method iteration?

here is a shooting method solution right out of the docs:

sol = First[
NDSolve[{x''[t] + Sin[x[t]] == 0 , x[0] == x[10] == 0}, x, t,
Method -> {"Shooting",
"StartingInitialConditions" -> { x'[0] == 1.666 }}]]
Plot[Evaluate[x[t] /. sol], {t, 0, 10}]


The shooting method is of course iterative, so how to monitor progress (initial condidion vs end condition )? I came up with this sort of hack approach and I wonder if there is a better way:

define a function that always returns zero, but saves what we want as a side effect:

zero[t_?NumericQ, x_, xp_] :=
(If[t == 0, xp0 = xp, If[t == 10, Sow[{xp0,x}]]]; 0)


now add as a term in the equation (no effect on solution since its always 0):

r = Reap[NDSolve[{x''[t] + Sin[x[t]] == zero[t, x[t], x'[t]] ,
x[0] == x[10] == 0}, x, t,
Method -> {"Shooting",
"StartingInitialConditions" -> { x'[0] == 1.666}}]][[2, 1]];

ListPlot[r, Frame -> True, FrameLabel -> {"x'[0]", "x[10]"},
PlotRange -> {{1.5, 2}, Automatic}]


is there a better way? I tried working with WhenEvent with no luck..

• I think ParametricNDSolve[] is way more convenient to use if you're performing manual shooting; that works with WhenEvent[] as well. Commented Jun 21, 2016 at 17:38
• I'm sure it is, but in this case I was trying to examine what the built in shooting method is doing. Commented Jun 21, 2016 at 17:42

Maybe commandeer FindRoot with the Villegas-Gayley trick: Updated, with the order of the steps taken by FindRoot saved in icsteps. The results of FindRoot, as saved by NDSolve and shown below as DownValues[], have been sorted by Mathematica and are not in the order in which there were called. This update stores the order in icsteps.

Clear[x, t];
InternalInheritedBlock[{FindRoot},
Unprotect[FindRoot];
FindRoot[f_, vars_, rest__] /; ! TrueQ[$in] := Block[{$in = True, res},
(* saved data: nf, objective, call, icsteps *)
nf = f;                                   (* NumericalFunction (unused here) *)
objective = DownValues @@ {Head@f["FunctionExpression"]}; (* expression used by f/nf *)
call = Inactive[FindRoot][f, vars, rest]; (* intended call, minus the monitor below *)
{res, icsteps} = Reap[
Sow[vars[[1, 1]], "icsteps"];
FindRoot[f, vars, rest, StepMonitor :> (Sow[#, "icsteps"] &)],
"icsteps"];
res];
Protect[FindRoot];
sol = First[
NDSolve[{x''[t] + Sin[x[t]] == 0, x[0] == x[10] == 0}, x, t,
Method -> {"Shooting", "StartingInitialConditions" -> {x'[0] == 1.666}}]];
]


The function passed to FindRoot:

objective


The actual call:

call
(*
Inactive[FindRoot][
ExperimentalNumericalFunction[{NDSolveShootingShootingDumpsic$762878}, NDSolveShootingShootingDumpg$762878[NDSolveShootingShootingDumpsic$762878], "-NumericalFunctionData-"], {{{0., 1.666}}}, "AccuracyGoal" -> 7.97729, "PrecisionGoal" -> 7.97729, "MaxIterations" -> Automatic, "Method" -> {"Newton", "StepControl" -> "TrustRegion"}, "WorkingPrecision" -> MachinePrecision] *)  Forgot to add the NDSolveReinitialize call and results, which are memoized: Cases[objective, HoldPattern[NDSolveShootingShootingDumpres$ = ndcall_Symbol[_Symbol]] :>
DownValues[ndcall], Infinity]
(*
{{HoldPattern[NDSolveShootingShootingDumpsol$767038[{-0.151661, 1.86146}]] :> {{{-0.151661, 1.86146}, {0.127858, 1.86324}}, {{{1., 0.}, {0., 1.}}, {{3.036, -25.0738}, {-0.288852, 2.71496}}}}, ... a bunch of results from trial integrations... HoldPattern[NDSolveShootingShootingDumpsol$767038[{0., 1.87817}]] :>
{{{0., 1.87817}, {-0.0000590907, 1.87817}},
{{{1., 0.}, {0., 1.}}, {{1., -27.9028}, {0.000031445, 0.999123}}}},
HoldPattern[
NDSolveShootingShootingDumpsol$767038[NDSolveShootingShootingDumpic$_?VectorQ]] :>
(NDSolveShootingShootingDumpsol$767038[NDSolveShootingShootingDumpic$] =
Module[{NDSolveShootingShootingDumptstate$, NDSolveShootingShootingDumpsres$},
NDSolveShootingShootingDumptstate$= Quiet[NDSolveReinitialize[ NDSolveShootingShootingDumpstate$767038,
{NDSolveShootingShootingDumpx$767038[NDSolveShootingShootingDumpts$767038] ==
NDSolveShootingShootingDumpic$, NDSolveShootingShootingDumpY$767038[NDSolveShootingShootingDumpts$767038] == NDSolveShootingShootingDumpid$767038}],
{NDSolveReinitialize::precw}];
If[! ListQ[NDSolveShootingShootingDumptstate$], Throw[$Failed]];
NDSolveShootingShootingDumptstate$= First[NDSolveShootingShootingDumptstate$];
NDSolveShootingShootingDumpsres$= NDSolveShootingShootingDumpiterate[NDSolveShootingShootingDumptstate$];
{NDSolveShootingShootingDumpsres$[[ All, 1 ;; NDSolveShootingShootingDumpn$767038
]],
(Partition[#1, NDSolveShootingShootingDumpn$767038] &) /@ NDSolveShootingShootingDumpsres$[[
All, NDSolveShootingShootingDumpn\$767038 + 1 ;; -1
]]
}])}}
*)


As noted above, these downvalues do not show the order of the steps. The steps may be seen in icsteps:

icsteps
(*
{{{0., 1.666}, {-0.00021824, 1.8383}, {-0.00124324, 1.83842},
{-0.00327712, 1.83866}, {-0.00727959, 1.83913}, {-0.0150193, 1.84003},
{-0.0294259, 1.84172}, {-0.0540885, 1.84465}, {-0.0974275, 1.8502},
{-0.151661, 1.86146}, {-0.0495448, 1.87856}, {-6.93889*10^-18, 1.87882},
{0., 1.87817}, {0., 1.87817}, {0., 1.87817}}}
*)

• Note: I orginally used Print instead of storing the arguments in variables. FindRoot is called only once. This is not clear from the present code. Commented Jun 21, 2016 at 18:24