3
$\begingroup$

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:

error message

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?

$\endgroup$
5
$\begingroup$

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.

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

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

Internal`AddHandler["Message", messageHandler]
(* lots of NDSolve executions *)
Internal`RemoveHandler["Message", messageHandler]

Or, if you want it more contained

Internal`HandlerBlock[{"Message", messageHandler},
  Message[NDSolve::ndsz, "T", 5]
]
vals
(* 5 *)
$\endgroup$
2
$\begingroup$

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"]]

enter image description here

$\endgroup$

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.