# Catching NDSolve Singularities

I am using NDSolve to solve a system of equations with changing initial conditions (I'm using NDSolve multiple times). However, for particular initial conditions I am getting an effective step size of 0 during the integration. Is there a way to detect this in NDSolve and output some useful info?

i.e.:

x10 = 1.2
x20 = 3.2
eqn1 = x1''[t] ==  (x1[t] - x2[t]) / (x1[t] - x2[t])^2
eqn2 = x2''[t] ==  (x1[t] - x2[t]) / (x1[t] - x2[t])^2

NDSolve[{eqn1, eqn2, x1[0]==x10, x2[0]==x20, x1'[0]==0, x2'[0]==0}, {x1[t], x2[t]}, {t,0,1}]


I have tried adding the following to NDSolve:

WhenEvent[(x1[t]-x2[t])==0, "StopIntegration"; Print[x10,x20]]


But the WhenEvent command doesn't stop the integration before it sees the singularity, and thus I don't see the Print command. I need to see the Print command because I have code that loops over the first block above and tries random values of x10 and x20. I need to know what initial conditions make the integration fail.

• Have you tried something like Check[NDSolve[..], {x10,x20}]? Jan 20, 2020 at 21:39

You can monitor the step size. It will be "effectively zero" well before it is zero. The code below does not quite do what you want, I think, because it cannot trap NDSolve after it has decided to exit. When NDSolve fails due to step size "zero," it does not take the last step. At best you can only catch the next-to-last step, and probably you have to aim a little earlier than that. But maybe it's helpful.

The threshold to use is one of the weak points. It depends on the Method. The "Extrapolation" method seems really hard to trap, because it can fail even when the next-to-last step is large.

SeedRandom[4];
xcurrent =.;
threshold = 75 + 1; (* relative size to x * $$MachineEpsilon *) y0 = 1.; sol = NDSolve[{y'[x] == y[x]^2, y[0] == y0, WhenEvent[ x - threshold*$$MachineEpsilon*Abs[x] - xcurrent < 0 // TrueQ,
{Print["Whoa, there!: ", {x, y0, y[x], NDSolveSelf}];
y[x] -> (y0 = RandomReal[])}]},
y, {x, 0, 10}, StepMonitor :> (xcurrent = x)]


Whoa, there!: {1.,1.,2.15996*10^12,NDSolveStateData[<0.>]}

Whoa, there!: {4.21155,0.311376,3.971*10^11,NDSolveStateData[<0.,1.>]}

Whoa, there!: {5.97832,0.566004,3.73592*10^11,NDSolveStateData[<0.,4.21155>]}

Whoa, there!: {7.57282,0.627155,2.85592*10^11,NDSolveStateData[<0.,5.97832>]}

Also beware the step size criterion is no good when x` is close to zero. One should replace it with an absolute step size appropriate to the problem at hand.