1

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?

Community
  • 1
dearN
  • 111
  • 2
  • 1
    Welcome to SciComp! Please include the ODE (and its initial conditions) using LaTeX syntax so that non-Mathematica users can better understand your example. – Geoff Oxberry Apr 04 '14 at 04:42
  • 1
    This might be helpful: http://blogs.mathworks.com/cleve/2014/06/09/ordinary-differential-equations-stiffness/?s_tid=Blog_Cleve_Category – Rhei Jan 20 '15 at 13:44
  • Someone (Gustaf Soderlind) has finally proposed a quantitative measure of stiffness that experts are beginning to agree on. See this paper. – David Ketcheson Jan 15 '17 at 08:07

1 Answers1

2

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.

Community
  • 1
Geoff Oxberry
  • 30,394
  • 9
  • 64
  • 127
  • Yes, I am currently trying to plot item 2 of your description. Will put that down here as an answer soon. – dearN Apr 04 '14 at 19:31