8
$\begingroup$

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

enter image description here

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

enter image description here

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

$\endgroup$
2
  • 1
    $\begingroup$ I think ParametricNDSolve[] is way more convenient to use if you're performing manual shooting; that works with WhenEvent[] as well. $\endgroup$ – J. M.'s torpor Jun 21 '16 at 17:38
  • $\begingroup$ I'm sure it is, but in this case I was trying to examine what the built in shooting method is doing. $\endgroup$ – george2079 Jun 21 '16 at 17:42
8
$\begingroup$

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

Mathematica graphics

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}}}
*)
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$ – Michael E2 Jun 21 '16 at 18:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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