inspecting step size and order of $\tt NDSolve$

Question

I am trying to collect information about what step sizes and what orders is using NDSolve internally. I tried wrapping it into a Trace but it dumping a huge amount of non-relevant info.

Any idea how to retrieve this easily?

Answer 1

There are a number of ways for each question.

The steps may be obtained from the solution itself.

ysol = First@
   NDSolve[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1}, y, {x, 0, 30}];

steps = y["Coordinates"] /. ysol // First // Differences;
steps // ListPlot

The use of MethodMonitor to extract method information is shown in user21's answer to How to find out which method Mathematica selected?. But while the example there,

Reap[sol = 
   NDSolve[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1}, y, {x, 0, 30}, 
    Method -> "StiffnessSwitching", 
    MethodMonitor :> (Sow[NDSolve`Self[[0]]];)];]

yields a lot of information about "StiffnessSwitching", applying it to the Automatic method yields none,

Reap[sol = 
   NDSolve[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1}, y, {x, 0, 30}, 
    MethodMonitor :> (Sow[NDSolve`Self[[0]]];)];]
(*  {Null, {}}  *)

One can use Trace. Here is a slightly difference pattern to search for than the one given in Simon's answer to the linked question.

Trace[
 NDSolve[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1}, y, {x, 0, 30}],
 NDSolve`InitializeMethod[meth_, data___] :> meth,
 TraceInternal -> True]
(*  {{{Automatic, NDSolve`LSODA}}}  *)

At first the method is initialized to Automatic, but when integration starts, it is set to LSODA.

Alternatively, one can use the tools for setting up analyzing NDSolve methods, which are described in the tutorials Components and Data Structures and NDSolve Method Plugin Framework.

You can set up an ODE with NDSolve`ProcessEquations, but the method won't be initialized yet:

{state} = NDSolve`ProcessEquations[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1}, 
   y, {x, 0, 30}];
state@"MethodData"["Forward"]
(*  None  *)

If you iterate (i.e., integrate over any time period), one can extract the method information. (Beware it only gives the current information.)

NDSolve`Iterate[state, 30]
state@"MethodData"["Forward"] // Short
(state@"MethodData"["Forward"])["DifferenceOrder"]
(*
  NDSolve`LSODA[NDSolve`MethodData[<<1>>]]
  8
*)

We see that the difference order is 8.

So user21's use of MethodMonitor is useful when the method switches, but one of the other approaches above should be used when MethodMonitor doesn't return any information. On the other hand, to get the "DifferenceOrder", I seem to have to get the method object from the NDSolve`StateData state; it's not clear to me whether that information may be obtained through MethodMonitor.

Answer 2

I think for the order, it depends on the method used in NDSolve. But for the step size we can get easily from StepMonitor . For example:

We can get all the point where a step is taken in the solution process

data = Reap[NDSolve[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1}, y, {x, 0, 30}, 
    StepMonitor :> Sow[{x, y[x]}]]][[2, 1]]

The step size can be plotted as

ListPlot[Differences[data][[All, 1]], AxesLabel -> {"step No.", "step size"}]