# Automatically detect largest interval over which NDSolve can find a solution

Question: Consider the following numerical resolution:

NDSolve[eqn, {x1[t], x2[t], y[t]}, {t, tmin, tmax}


where eqn, tmin and tmax are given.

Suppose this gives an error (for example, a stiff system detection, or any numerical error).

How would you do to adjust automatically tmin and tmax to the largest (or a "large") interval where NDSolve doesn't give any error?

I tried using Check but I don't manage to use it efficiently as I cannot catch the defective times.

Context: In case you're wondering, I'm asking because I have a large number of odes eqn to solve, some of which work fine on $[t_\mathrm{min},t_\mathrm{max}]$, but some others don't. And of course I don't want to adjust the time interval by hand!

Example: With

eqn = {x1'[t] == -x1[t]^2 - x2[t] + y[t]^3, x2'[t] == x1[t] - x2[t] + x1[t]^2 x2[t]^2,
y'[t] == x2[t]^2 + y[t] + x1[t]^2 y[t]^2 + 0.5 y[t]^3, x1[0] == 1, x2[0] == 1, y[0] == 1};


the following

{tmin, tmax} = {-1, 1};
NDSolve[eqn, {x1[t], x2[t], y[t]}, {t, tmin, tmax}];


yield an NDSolve::ndsz error while

{tmin, tmax} = {-0.5, 0.5};
NDSolve[eqn, {x1[t], x2[t], y[t]}, {t, tmin, tmax}];


works fine.

• Actually, eqn becomes stiff at t~0.235. Have you tried Eventlocator? – Feyre Aug 10 '14 at 9:03
• I have, i did not managed to make it general (is it an upper or a lower bound, etc.). – anderstood Aug 10 '14 at 16:34

NDSolve has already detected the largest such intervals for you, which is why the resulting InterpolatingFunctions have restricted domains. You can use InterpolatingFunctionDomain to extract those domains. I'd do something like so

Clear[x1, x2, y]
eqn = {x1'[t] == -x1[t]^2 - x2[t] + y[t]^3,
x2'[t] == x1[t] - x2[t] + x1[t]^2 x2[t]^2,
y'[t] == x2[t]^2 + y[t] + x1[t]^2 y[t]^2 + 0.5 y[t]^3, x1[0] == 1,
x2[0] == 1, y[0] == 1};
Quiet[
{x1, x2, y} = NDSolveValue[eqn, {x1, x2, y}, {t, -1, 1}],
NDSolveValue::ndsz]


Note that the domains are indicated as part of the result. If you need to automatically extract them within code, you can do something like so

Needs["DifferentialEquationsInterpolatingFunctionAnatomy"];
domains = InterpolatingFunctionDomain /@ {x1, x2, y}
{a, b} = {Min[#], Max[#]} &[Intersection[Interval @@@ domains]]
(* Out:
{{{-0.720382, 0.234762}}, {{-0.720382, 0.234762}}, {{-0.720382, 0.234762}}}

{-0.720382, 0.234762}
*)


Of course, the point is that you can use this in subsequent command without worrying about the problems arising from trying to use inputs from outside the domain of the InterpolatingFunctions. For example, you can generate a plot.

ParametricPlot[{{x1[t], y[t]}, {x2[t], y[t]}}, {t, a, b}]


• +1. In a single NDSolve call like this, can't one assume that the domains of x1, x2, y will all be the same? – Michael E2 Aug 10 '14 at 14:52
• Perfect, with a few adjustements it worked (I'm on version 7), thank you. @MichaelE2: if think we can, and we save one line :) – anderstood Aug 10 '14 at 16:38
• @MichaelE2 You are probably right - I was clearly being paranoid. :) – Mark McClure Aug 11 '14 at 1:33

This seems to work for me:

thisstep = 0;
laststep = 0;
eqn = {x1'[t] == -x1[t]^2 - x2[t] + y[t]^3,
x2'[t] == x1[t] - x2[t] + x1[t]^2 x2[t]^2,
y'[t] == x2[t]^2 + y[t] + x1[t]^2 y[t]^2 + 0.5 y[t]^3, x1[0] == 1,
x2[0] == 1, y[0] == 1};
{tmin, tmax} = {-1,1};

First@NDSolve[eqn, {x1[t], x2[t], y[t]}, {t, tmin, tmax},
MaxStepFraction -> 1/150,
StepMonitor :> (laststep = thisstep; thisstep = t;
stepsize = thisstep - laststep;),
Method -> {"MethodOfLines",
Method -> {"EventLocator",
"Event" :> (If[stepsize < 10^-4, 0, 1])}}
];


You might have to play around with minimum stepsize.

• It yields the error NDSolve::molpde` (because of "MethodOfLines"). But anyway I think it is easier to simply take the bounds from the output of NDSolve. Since this works, and do not investigate further. Thank you for your suggestion. – anderstood Aug 10 '14 at 16:42