Does NDSolve iterate faster in a region where the system being solved is in equilibrium?

I was wondering how exactly the time it takes for NDSolve to iterate through a certain region of the domain depends on the local behaviour of the given system being solved. My guess would be the greater in absolute value the derivatives of the system are in the region (with respect to the domain variable), the smaller the steps NDSolve takes, and thus the more computing time it needs, but I'm not sure, and I'm particularly unsure about the time it takes for NDSolve to iterate through a region where the system is in equilibrium. Is it costly?

• It's a rather broad question. If an adaptive-step method is used, the step size is adjusted so that the local error estimate is sufficiently small. The steps thus will depend on the order & error estimate of the method. If the error is very small, NDSolve will increase the step size every few steps (by up to a factor of 10, apparently). – Michael E2 May 1 '16 at 22:13

Here is a simple example for the default method (LSODA) that shows a couple of the issues related to the situation in the question:

$$y' = \exp(-10 \, x^2), \quad y(-20) = 1, \quad -20 \le x \le 20$$

The exact solution is a scaled error function (Erf), which has a sigmoid graph going from $y(-20)=1$ to $y(20) \approx 1.5605$. For $|x|$ large, $y'$ is extremely close to zero. I'll assume that such regions will serve as regions where the system is in "equilibrium," even though it is approximate.

The command ListPlot[sol] of an InterpolatingFunction solution sol shows the steps NDSolve takes. One can get the numerical values of the steps with sol["Grid"]. There are options MaxStepSize and MaxStepFraction that limit how large a step may be taken. By default MaxStepFraction is 1/10. These can be adjusted to suit the scale of an ODE.

The default behavior on our example is disappointing. NDSolve has been too fast and stepped right over the active region from one equilibrium region to the other.

x0 = 20;
yFinal = 1 + Integrate[Exp[-10 x^2], {x, -x0, x0}];
sol = NDSolveValue[{y'[x] == Exp[-10 x^2], y[-x0] == 1}, y, {x, -x0, x0}];
ListPlot[
sol,
GridLines -> {None, {yFinal}},
PlotRange -> {{-x0, x0}, {1 - 0.02, yFinal + 0.02}}]


We'll come back to this, but reducing the step size cures the problem:

sol2 = NDSolveValue[{y'[x] == Exp[-10 x^2], y[-x0] == 1}, y, {x, -x0, x0},
MaxStepFraction -> 1/20];
ListPlot[
sol2,
GridLines -> {None, {yFinal}},
PlotRange -> {{-x0, x0}, {1 - 0.02, yFinal + 0.02}}]


There are couple of things to note about the plots. The maximum steps sizes are 4, 2 (resp.), corresponding to MaxStepFraction settings of 1/10, 1/20 (resp.) of an interval [-20,20] of length 40. We can see the maximum step being achieved in the equilibrium regions. It's not clear from the graph, but the steps do not start out at the maximum, which we will show. There is also a curious change of step size around x = -4 in the first plot, which we will explain.

Let's examine the result of the first code. Here are the actual steps:

sol["Grid"]
(*
{{-20.}, {-19.9979}, {-19.9959}, {-15.9959}, {-11.9959}, {-7.99589},
{-3.99589}, {-3.99555}, {-3.99521}, {-3.99453}, {-3.99384},
{-3.99316}, {-3.98631}, {-3.97947}, {-3.97263}, {-3.90419},
{-3.83575}, {-3.76731}, {-3.08293}, {-2.39855}, {-1.71417},
{2.28583}, {6.28583}, {10.2858}, {14.2858}, {17.1429}, {20.}}
*)


The first few steps are usually quite small. These are used to get an approximation of the local error, and then the real initial step size is computed. The function is so flat that the maximum step size 4. is selected. We can see the step size more clearly with this:

Differences@sol["Grid"]
(*
{{0.00205297}, {0.00205297}, {4.}, {4.}, {4.}, {4.}, {0.00034219},
{0.00034219}, {0.000684381}, {0.000684381}, {0.000684381},
{0.00684381}, {0.00684381}, {0.00684381}, {0.0684381}, {0.0684381},
{0.0684381}, {0.684381}, {0.684381}, {0.684381}, {4.}, {4.}, {4.},
{4.}, {2.85708}, {2.85708}}
*)


Now what happens at the seventh step x == -3.99589? The step size suddenly gets much smaller -- why? The next step would be at this number plus 4., but if we evaluate the right-hand side, we (or at least NDSolve) gets a surprise:

Exp[-10 x^2] /. x -> sol["Grid"][[7]] + 4.
(*  {0.999831}  *)


Suddenly the RHS is quite large by comparison. So the step is rejected, and NDSolve restarts the integration with a very small step. At this point, NDSolve is somewhat cautious, increasing the step size every few steps, first by a factor of 2 and then by a factor of 10, until the maximum step size 4. is reached again. At this point, it is a bit closer to the next equilibrium domain:

sol["Grid"][[21]]
(*  {-1.71417}  *)


The next step size is 4. and at that point, the RHS evaluates to a small number:

Exp[-10 x^2] /. x -> sol["Grid"][[21]] + 4.
(*  {2.03232*10^-23}  *)


As a result, this step is accepted and the region in the middle where the RHS is large is skipped over; NDSolve continues until it gets near the end. (The last couple steps normally subdivide the remainder of the interval.)

In the second code, the integration is restarted near x == -2 for the same reason as the first case; however, the initial steps from that point carry the integration into a region where the RHS is relatively large. Consequently, the maximum step size is not reached again until the next equilibrium region is reached.

• that was something I was not expecting to learn and totally worth reading. Thanks for the detailed explanation. – tsuresuregusa May 2 '16 at 20:13