# Is it possible to assign numbers generated in an Error-Message?

If I have a differential equation of the form:

NDSolve[{...}, {x, y, z}, {t, 0, 2}, WorkingPrecision -> 30,
MaxSteps -> Infinity, Method -> Automatic, InterpolationOrder -> All]


and I get a singularity at a specific point of t, for example if an object falls into a black hole, which leads to an error message like this: Is it possible to assign that value, here 1.2574' to a variable automatically, so that I can plot exactly until this point is reached?

If I just copy and paste that value manually I often end up a bit too early or too late, maybe because of roundig errors. Is there a way to automatically assign some variable χ to it, so that I can do a

Plot[..., {t, 0, χ}]


or do I have to copy and paste, and then play around with the last digits of the error message until I get as close to the singularity as possible?

Albert Retey has demonstrated in a similar situation that you can use "EventLocator" to detect an event in NDSolve. For example:

eqn = {\!$$\*SubscriptBox[\(∂$$, $$t$$]$$u[t, x]$$\) == 1/100 \!$$\*SubscriptBox[\(∂$$, $$x, x$$]$$u[t, x]$$\) - u[t, x] \!$$\*SubscriptBox[\(∂$$, $$x$$]$$u[t, x]$$\),
u[0, x] == Sin[2 π x], u[t, 0] == u[t, 1]};

NDSolve[eqn, u, {t, 0, 2}, {x, 0, 1}]


NDSolve::eerr: Warning: scaled local spatial error estimate of 5.741306825597143*^13 at t = 0.4450518534682055 in the direction of independent variable x is much greater than the prescribed error tolerance.

When the stiffness happens, Mathematica will try to take an effectively zero stepsize. You can see that by

NDSolve[eqn, u, {t, 0, 2}, {x, 0, 1},
StepMonitor :> (laststep = thisstep; thisstep = t;
stepsize = thisstep - laststep; Print[stepsize];)]
(*
0.0000115314
0.0000115314
8.70237*10^-6
...
...
7.88258*10^-15
*)


So we can use the small step size as a criteria to test the stiffness, and stop the integration

 NDSolve[eqn, u, {t, 0, 2}, {x, 0, 1},
StepMonitor :> (laststep = thisstep; thisstep = t;stepsize = thisstep - laststep;),
Method -> {"EventLocator", "Event" :> (If[stepsize < 10^-4, 0, 1])}]


Then the integration will stop when the step size is less than 10^-4, and the variable thisstep will be the point you are looking for.

• wow, using StepMonitor for its intended purpose. I wish I had thought of that. +1 – rcollyer May 27 '16 at 2:19

A method is to write a message handler, like in this answer. The handler is passed an argument of the form

Hold[Message[...], boolean]


where the boolean tells the handler if the message is to be displayed, or not. Since you are looking to capture the info passed to NDSolve::ndsz, I would write the handler like

Clear[messageHandler, vals];
vals = {};
messageHandler[Hold[Message[NDSolve::ndsz, _, v_], True]] :=
(vals = Flatten[{vals, v}])


Then, if you have multiple executions of NDSolve to perform, I would use

InternalAddHandler["Message", messageHandler]
(* lots of NDSolve executions *)
InternalRemoveHandler["Message", messageHandler]


Or, if you want it more contained

InternalHandlerBlock[{"Message", messageHandler},
Message[NDSolve::ndsz, "T", 5]
]
vals
(* 5 *)


Why not just extract the domain from the interpolating function? Using @xslittlegrass' example:

eqn = {
D[u[t, x], t] == 1/100 D[u[t,x], x, x] - u[t, x] D[u[t, x], x],
u[0, x] == Sin[2 π x], u[t, 0] == u[t, 1]
};

sol = NDSolveValue[eqn, u, {t, 0, 2}, {x, 0, 1}];


NDSolveValue::ndsz: At t == 0.44505185346820575, step size is effectively zero; singularity or stiff system suspected.

NDSolveValue::eerr: Warning: scaled local spatial error estimate of 5.741306825596747*^13 at t = 0.44505185346820575 in the direction of independent variable x is much greater than the prescribed error tolerance. Grid spacing with 25 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or a smaller grid spacing can be specified using the MaxStepSize or MinPoints method options.

The domain is then:

sol["Domain"]


{{0., 0.445052}, {0., 1.}}

and the desired value is:

sol["Domain"][[1,2]]


0.445052

Alternatively, you can construct a region from the domain, and use this in your plotting function:

Plot3D[sol[x, t], {x, t} ∈ Rectangle @@ Transpose @ sol["Domain"]] 