2
$\begingroup$

I want to detect directly when step size becomes "effectively" zero. The below example from the documentation throws the error message as expected:

s = {};
NDSolve[{(2 - f[x]) f'[x] == f[x], f[0] == 1}, f, {x, 0, 5}, StepMonitor :> AppendTo[s, x]];

NDSolve::ndsz: At x == 0.3862940268757776`, step size is effectively zero; singularity or stiff system suspected.

The below code indicates that none of the actual steps taken has zero length.

AnyTrue[Differences@s, PossibleZeroQ]

(* False *)

How does NDSolve decide that the step size is zero? I can of course capture the NDSolveValue::ndsz error, but I want to know when exactly (depending on what parameters) the error is issued. In some extreme cases, NDSolve may generate an InterpolatingFunction solution that has practically zero domain length (but not according to PossibleZeroQ).

$\endgroup$
1
  • 1
    $\begingroup$ I think step size is taken to be zero, if the difference between two steps is nearly working precision. Look at: Differences@ Reap[NDSolve[{(2 - f[x]) f'[x] == f[x], f[0] == 1}, f, {x, 0, 5}, StepMonitor :> Sow[{x, f[x]}], WorkingPrecision -> 30]][[2, 1, -10 ;; -1]] and the same without WorkingPrecision->... $\endgroup$ – Akku14 Oct 13 '20 at 15:29
3
$\begingroup$

Update 2021.01.09: I figured out how to examine the step size that is effectively zero.

The NDSolve::ndsz results when adding the next step h to the current value of x results in a value equal to x. Since Equal is compared with tolerance, adding h is like adding zero. Hence the "step size is effectively zero." (This tolerance does not depend on Internal`$EqualTolerance. Either NDSolve resets Internal`$EqualTolerance or it is internally hard-coded. I've tested several examples, and the behavior is consistent with a tolerance equal to the default value of Internal`$EqualTolerance.)

{state} = 
  NDSolve`ProcessEquations[{(2 - f[x]) f'[x] == f[x], f[0] == 1}, 
   f, {x, 0, 5}];
NDSolve`Iterate[state, 5];
solIF = f /. NDSolve`ProcessSolutions[state];

NDSolve`Iterate::ndsz : At x == 0.3862940268757776, step size is effectively zero; singularity or stiff system suspected.

lastcoord = solIF@"Coordinates" // First // Last
nexth = state@"TimeStep"["Forward"]
lastcoord == lastcoord + nexth
(*
  0.386294        <-- x
  3.35604*10^-15  <-- h
  True            <-- x == x + h
*)

Original answer:

A short answer to support Akku14's remark: "...a step size that is too small for numerically reliable calculations." (docs for NDSolve::ndsz).

The last step when NDSolve stops (in this case, with LSODA method) is usually a few hundred times the relative epsilon for the working precision ($\approx 2 \times 10^{-p}$). I've seen a wide range from less than ten to almost 1000 times epsilon.

sol = NDSolve[{(2 - f[x]) f'[x] == f[x], f[0] == 1}, f, {x, 0, 5}];
With[{steps = f["Grid"] /. sol // Flatten},
 Last@Differences[steps]/(Last@steps*2*10^-Precision[sol])]
(*  385.069  *)

I think, at least with LSODA, the problem is that the following is probably true:

You do not see the step size that causes the integration to stop.

It's the next step that would be effectively zero. One can see the stiffness develop in the rather rapid decrease of step size, the large gaps being from repeated error test failures:

With[{steps = f["Grid"] /. sol // Flatten},
  Differences[steps]] // ListLogPlot

I don't know how (or if) you can get out of LSODA data about its current state. You could test when the steps size gets below 1000 times epsilon. That seems to be a possible threshold.

$\endgroup$
3
  • $\begingroup$ This all sounds very reasonable, thanks for the information. I will run some tests on my use case. To be honest, I am a bit disappointed running into some magic numbering again. Unreliable factors like these decrease robustness of code. $\endgroup$ – István Zachar Oct 13 '20 at 21:00
  • $\begingroup$ @IstvánZachar If I read the ODEPACK code correctly, the step size in LSODA might be reduced by a factor as small as $10^{-8}$ without generating a NDSolve::nderr or NDSolve::ndcf error. Conceivably the "next step" that is too small might be on the order of epsilon or less; also "1000 times epsilon" as a threshold should fail to prevent NDSolve::ndsz sometimes. $\endgroup$ – Michael E2 Oct 14 '20 at 18:30
  • $\begingroup$ Thanks Michael, this is a nice explanation now. $\endgroup$ – István Zachar Jan 9 at 13:53

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.