I am using NDSolver and WhenEvent to solve a system of differential equations. I get this message:

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

The solution I get is pretty stable until that time, then blows off completely. At each evaluation I get the same behaviour at different times.

I understand that the problem might lie with the fact that I am solving for very long times, and eventually the algorithm just starts having problems handling all the data.

Is there a way to stop the solver as soon as the error comes up?

  • 1
    $\begingroup$ Doesn't the solver already stop exactly at the point mentioned? $\endgroup$ – J. M.'s discontentment Jul 4 '15 at 10:52
  • $\begingroup$ Maybe it does, but the solution is then just extrapolated for the remaining time. Is there a way to return the last time t before the error? I am thinking of something like WhenEvent["error", tfinal = t] $\endgroup$ – Andrea Jul 4 '15 at 11:29
  • $\begingroup$ Maybe this mathematica.stackexchange.com/questions/66590/… $\endgroup$ – Mariusz Iwaniuk Jul 4 '15 at 13:13
  • $\begingroup$ Unfortunately the system does not hit a singularity, so I have no way to analytically identify a point where to stop the integration, like in the link. There is no way to make it stop in response to the error itself? $\endgroup$ – Andrea Jul 4 '15 at 13:28
  • 3
    $\begingroup$ Let's say you have sol = NDSolve[..., x, {t, 0, 100000}]. The x["Domain"] /. sol will yield {{{0, 49169.954923435114}}}, and you can get the time it stopped. You can also use the option NDSolve[..., "ExtrapolationHandler" -> {Indeterminate&}] and the solution won't extrapolate; however, it will return Indeterminate, which may be good or may cause other problems. Without explicit code, we're just guessing. $\endgroup$ – Michael E2 Jul 4 '15 at 17:43

First, a test case:

sol = NDSolve[{x'[t] == -1/x[t], x[0] == -1}, x, {t, 0, 2}]

NDSolve::ndsz: At t == 0.499999735845871`, step size is effectively zero; singularity or stiff system suspected. >>

Mathematica graphics

Consult What's inside InterpolatingFunction[{{1., 4.}}, <>]? or InterpolatingFunctionAnatomy and you will discover that InterpolatingFunction comes with lots of utilities and methods that give access to a wealth of information. In particular, the method "Domain" will yield the domain of the InterpolatingFunction, which in this case will be the interval over which NDSolve successfully integrated the ODE (up to the point where stiffness stopped it).

{tstart, tstop} = Flatten[x["Domain"] /. sol]
(*  {0.`, 0.499999735845871`}  *)

A couple of ways to plot. Note for ListLinePlot, you don't need to know the domain; it and ListPlot handle InterpolatingFunction as a special case.

ListLinePlot[x /. sol]
Plot @@ {{x[t], x'[t]} /. First[sol], Flatten@{t, x["Domain"] /. sol}}

Mathematica graphics

By default, going past the ends of the domain results in extrapolation.

Plot[x[t] /. First[sol], {t, 0, 2}]

Mathematica graphics

One can use the option, mentioned in a comment in the linked question, "ExtrapolationHandler" to override the default extrapolation.

sol2 = NDSolve[{x'[t] == -1/x[t], x[0] == -1}, x, {t, 0, 10}, 
   "ExtrapolationHandler" -> {Indeterminate &}];

Plot[x[t] /. First[sol2], {t, 0, 2}]

NDSolve::ndsz: At t == 0.4999997164125506`....

Mathematica graphics

| improve this answer | |

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.