NDSolve how to monitor shooting method iteration?

Question

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..

Michael E2 · Answer

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];
Internal`InheritedBlock[{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][
 Experimental`NumericalFunction[{NDSolve`Shooting`ShootingDump`sic$762878},
  NDSolve`Shooting`ShootingDump`g$762878[NDSolve`Shooting`ShootingDump`sic$762878], 
  "-NumericalFunctionData-"], {{{0., 1.666}}}, 
 "AccuracyGoal" -> 7.97729, "PrecisionGoal" -> 7.97729, 
 "MaxIterations" -> Automatic, 
 "Method" -> {"Newton", "StepControl" -> "TrustRegion"}, 
 "WorkingPrecision" -> MachinePrecision]
*)

Forgot to add the NDSolve`Reinitialize call and results, which are memoized:

Cases[objective, 
 HoldPattern[NDSolve`Shooting`ShootingDump`res$ = ndcall_Symbol[_Symbol]] :> 
  DownValues[ndcall], Infinity]
(*
{{HoldPattern[NDSolve`Shooting`ShootingDump`sol$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[NDSolve`Shooting`ShootingDump`sol$767038[{0., 1.87817}]] :>
     {{{0., 1.87817}, {-0.0000590907, 1.87817}},
      {{{1., 0.}, {0., 1.}}, {{1., -27.9028}, {0.000031445, 0.999123}}}}, 
  HoldPattern[
    NDSolve`Shooting`ShootingDump`sol$767038[NDSolve`Shooting`ShootingDump`ic$_?VectorQ]] :> 
     (NDSolve`Shooting`ShootingDump`sol$767038[NDSolve`Shooting`ShootingDump`ic$] = 
     Module[{NDSolve`Shooting`ShootingDump`tstate$, 
       NDSolve`Shooting`ShootingDump`sres$}, 
      NDSolve`Shooting`ShootingDump`tstate$ = 
       Quiet[NDSolve`Reinitialize[
         NDSolve`Shooting`ShootingDump`state$767038,
         {NDSolve`Shooting`ShootingDump`x$767038[NDSolve`Shooting`ShootingDump`ts$767038] == 
           NDSolve`Shooting`ShootingDump`ic$, 
          NDSolve`Shooting`ShootingDump`Y$767038[NDSolve`Shooting`ShootingDump`ts$767038] == 
           NDSolve`Shooting`ShootingDump`id$767038}],
     {NDSolve`Reinitialize::precw}]; 
      If[! ListQ[NDSolve`Shooting`ShootingDump`tstate$], Throw[$Failed]]; 
      NDSolve`Shooting`ShootingDump`tstate$ = First[NDSolve`Shooting`ShootingDump`tstate$]; 
      NDSolve`Shooting`ShootingDump`sres$ = 
       NDSolve`Shooting`ShootingDump`iterate[NDSolve`Shooting`ShootingDump`tstate$];
      {NDSolve`Shooting`ShootingDump`sres$[[
        All, 1 ;; NDSolve`Shooting`ShootingDump`n$767038
        ]],
       (Partition[#1, NDSolve`Shooting`ShootingDump`n$767038] &) /@ 
        NDSolve`Shooting`ShootingDump`sres$[[
         All, NDSolve`Shooting`ShootingDump`n$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}}}
*)

NDSolve how to monitor shooting method iteration?

One Answer

Add your own answers!

Ask a Question