1
$\begingroup$

On a recommendation from Mathematica.SE, I am posting this on Computational Science.SE:

I am trying to quantify stiffness of an ODE by relating it to the fine-ness with which NDSolve treats it's time step.

BACKGROUND

For the uninitiated:

Essentially, when an ordinary differential equation is stiff and we attempt to solve it with NDSolve in Mathematica, it would employ a stiff solver to detect and reconcile stiffness.

I use the stiff equation shown on wikipedia to demonstrate this. This equation is solved using BDF (backward difference formulation/formula) which is a recognized method to solve stiff equations.

Equation and solution using NDSolve with BDF

tMax=100;
rSol=
r/.NDSolve[
r'[t]==-15 r[t]&&r[0]==1,r,
{t,0,tMax},Method->{"BDF","MaxDifferenceOrder"->5}][[1]]

Plots of solution r vs t

{Plot[
rSol[t],
{t,0,100},
PlotRange->{{0,tMax},Automatic}
],
LogPlot[
rSol[t],
{t,0,100},
PlotRange->{{0,tMax},Automatic}
]}

Fig 1

enter image description here

Fig 2

enter image description here

Fig 1 shows that up to t=29.9 or so, there is considerable stiffness (unstable oscillations) that is reconciled by ensuring that the time step is sufficiently small.

THE REAL QUESTION:

What would be a good plot to quantify this stiffness and internal reconciliation/changing of time step? How can I best quantify stiffness graphically in Mathematica?

$\endgroup$
  • 1
    $\begingroup$ Welcome to SciComp! Please include the ODE (and its initial conditions) using LaTeX syntax so that non-Mathematica users can better understand your example. $\endgroup$ – Geoff Oxberry Apr 4 '14 at 4:42
  • 1
    $\begingroup$ This might be helpful: blogs.mathworks.com/cleve/2014/06/09/… $\endgroup$ – Rhei Jan 20 '15 at 13:44
  • $\begingroup$ Someone (Gustaf Soderlind) has finally proposed a quantitative measure of stiffness that experts are beginning to agree on. See this paper. $\endgroup$ – David Ketcheson Jan 15 '17 at 8:07
1
$\begingroup$

Stiffness has no precise quantitative definition. (See my answer here.)

Common numerical proxies would be things like:

  • stiffness ratio: find the (unsigned) ratio of largest magnitude eigenvalue to smallest magnitude eigenvalue of the Jacobian matrix of your right-hand side; estimates are probably fine
  • ratio of time step taken to maximum stable time step: a ratio closer to 1 indicates a stability-limited time step, so your problem is probably stiff; a ratio closer to 0 indicates an accuracy-limited time step, and probably a non-stiff problem

There may be other metrics out there; these are the first two that come to mind.

$\endgroup$
  • $\begingroup$ Yes, I am currently trying to plot item 2 of your description. Will put that down here as an answer soon. $\endgroup$ – dearN Apr 4 '14 at 19:31

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.