How to get the derivative with StepMonitor?

Question

I am trying to use StepMonitor inside of the NDSolveValue function for a complicated system of differential equations, but the simple example below shows my issue just fine.

I am trying to obtain values of the derivative during the process of solving the differential equation by using StepMonitor.

Reap[NDSolveValue[{y'[x] == -x, y[0] == 1}, y[x], {x, 0, 100},StepMonitor :> Sow[y'[x]]]]//Last

And I get as an output a list of values of derivatives at the various time points such as $y'[.00017]$ rather than what the derivative actually evaluates to at those time steps.

How do I get the actual value of the derivative rather than these expressions? Note that this needs to be done while NDSolveValue runs. I know I can get these values after the fact using the solution to the differential equation.

Answer 1

One possibility is to increase the order of your ODE:

Reap[
    NDSolveValue[
        {
        y''[x] == -1, (* obtained by differentiating y'[x] == -x *)
        y[0] == 1,
        y'[0] == 0
        },
        y[x],
        {x, 0, 100},
        StepMonitor :> Sow[y'[x]]
    ]
] //Last

{{-0.000145152, -0.000290305, -1.45181, -2.90334, -4.35486, -14.3549, -24.3549, -34.3549, -44.3549, -54.3549, -64.3549, -74.3549, -84.3549, -92.1774, -100.}}

Another possibility is to include an additional dependent variable (in which case NDSolveValue uses a DAE solver):

Reap[
    NDSolveValue[
        {
        y'[x] == z[x] == -x,
        y[0] == 1
        },
        y,
        {x, 0, 100},
        StepMonitor :> Sow[y'[x]]
    ]
] //Last

{{-0.0001, -0.0002, -0.00040102, -0.000711913, -0.001024, -0.00164782, -0.00227173, -0.00351954, -0.00601518, -0.0110064, -0.020989, -0.040954, -0.0808842, -0.160744, -0.320465, -0.639906, -0.959347, -1.59823, -2.87599, -5.43152, -10.5426, -20.5426, -30.5426, -40.5426, -50.5426, -60.5426, -70.5426, -80.5426, -90.5426, -100.}}

Answer 2

Breaking down NDSolveValue into its component steps and inserting a hook back to the numerical function that computes the derivative:

({state} = 
   NDSolve`ProcessEquations[{y'[x] == -x, y[0] == 1}, 
    y[x], {x, 0, 100}, StepMonitor :> Sow[First@nf[x, {y[x]}]]];
  nf = state@"NumericalFunction";
  NDSolve`Iterate[state, 100];
  y[x] /. NDSolve`ProcessSolutions[state]) // Reap
(*
  InterpolatingFunction[{{0., 100.}}, <<4>>][x],

  {{-0.000171081, -0.000342162, -0.000684323, -0.00102648, \
    -0.00136865, -0.00479026, -0.00821188, -0.0116335, -0.0458497, \
    -0.0800658, -0.114282, -0.456444, -0.798605, -1.14077, -4.56238, \
    -7.984, -11.4056, -21.4056, -31.4056, -41.4056, -51.4056, -61.4056, \
    -71.4056, -81.4056, -90.7028, -100.}}}
*)

---

Remark: First@nf[x, {y[x]}] returns the first element of the derivative vector, which is of the form {y'[x]} in this case. For higher order equations and higher dimensional systems, omitting First will give the derivatives of all the variables. Here's an example that uses the built-in utilities described in the tutorial Components and Data Structures to construct a step monitor that does not depend on how NDSolve orders the values in the derivative vector. It returns substitution rules that may be used to get whichever derivative values are desired.

({state} = 
   NDSolve`ProcessEquations[{y''[x] == -x, y[0] == 1, y'[0] == 0}, 
    y[x], {x, 0, 100},
    StepMonitor :> Sow[Thread[
       NDSolve`SolutionDataComponent[state@"Variables", "X'"] ->
        nf[
         NDSolve`SolutionDataComponent[state@"Variables", "T"],
         Through[
          NDSolve`SolutionDataComponent[state@"Variables", "X"]@
           NDSolve`SolutionDataComponent[state@"Variables", "T"]
          ]
         ]]]];
  nf = state@"NumericalFunction";
  NDSolve`Iterate[state, 100];
  y[x] /. NDSolve`ProcessSolutions[state]) // Reap // Last
(* omit "// Last" to get solution and derivative values
  {{{y' -> -2.06958*10^-8, y''] -> -0.00014386},
    {y' -> -6.20874*10^-8, y'' -> -0.000287721},
    ...
    {y' -> -5000., y'' -> -100.}}}
*)

Or to get just y''[x] and z''[x], in that order, from a system of second-order equations:

Clear[state, nf, ypp, zpp];
({state} = 
   NDSolve`ProcessEquations[{y''[x] == -x, y[0] == 1, y'[0] == 0, 
     z''[x] == -y[x] - 10, z[0] == 0, z'[0] == 0},
    {y[x], z[x]}, {x, 0, 100},
    StepMonitor :> Sow[Extract[  (* extract second derivatives from nf[] *)
       nf[
        NDSolve`SolutionDataComponent[state@"Variables", "T"],
        Through[
         NDSolve`SolutionDataComponent[state@"Variables", "X"]@
          NDSolve`SolutionDataComponent[state@"Variables", "T"]
         ]],
       Join[ypp, zpp]]]];
  ypp = Position[NDSolve`SolutionDataComponent[state@"Variables", "X'"], y''];
  zpp = Position[NDSolve`SolutionDataComponent[state@"Variables", "X'"], z''];
  nf = state@"NumericalFunction";
  NDSolve`Iterate[state, 100];
  {y[x], z[x]} /. NDSolve`ProcessSolutions[state]) // Reap