0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ You might also consider the MonitorMethod of the plug-in framework $\endgroup$ – Michael E2 Apr 7 at 15:19
  • $\begingroup$ Any feedback on my second answer? $\endgroup$ – Michael E2 Jun 21 at 14:26
  • $\begingroup$ @MichaelE2 I ended up answering the question I was trying to answer with this technique by taking a different approach, so I have to admit this fell to the wayside (it took me a bit to remember why I made this question in the first place). However, I do recognize I shouldn't leave questions/those who answer my questions hanging without feedback. When I have the time I will review your second answer in my system. $\endgroup$ – Aaron Stevens Jun 21 at 14:33
  • $\begingroup$ You ask a question, you don't explain a topic. $\endgroup$ – user64494 Jun 21 at 17:01
  • $\begingroup$ @user64494 What? $\endgroup$ – Aaron Stevens Jun 21 at 17:23
1
$\begingroup$

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

$\endgroup$
  • $\begingroup$ Adding more variables/equations will be harder/unreasonable for my actual application of this unfortunately. $\endgroup$ – Aaron Stevens Apr 4 at 20:46
1
$\begingroup$

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
$\endgroup$

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.