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$
3
  • 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$ Apr 4, 2014 at 4:42
  • 1
    $\begingroup$ This might be helpful: blogs.mathworks.com/cleve/2014/06/09/… $\endgroup$
    – Rhei
    Jan 20, 2015 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$ Jan 15, 2017 at 8:07

1 Answer 1

2
$\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$
1
  • $\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, 2014 at 19:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.