Can I use an explicit time stepping scheme to determine numerically whether an ODE is stiff?

Question

I have an ODE:

$u'=-1000u+sin(t)$
$u(0)=-\frac{1}{1000001}$

I know that this particular ODE is stiff, analytically. I also know that if we use an explicit (forward) time stepping method (Euler, Runge-Kutta, Adams, etc.), the method should return very large errors if the time step is too large. So, I have two questions:

1. Is this how stiff ODEs are determined, in general, when an analytical expression for the error term is not available or derivable?

2. In general, when the ODE is stiff, how do I determine a "small enough" timestep?

Answer 1

To answer your questions:

1. As far as I know, in practice, if explicit methods require extraordinarily small time steps relative to your time scale of interest (see answers to this question on what it means for an ODE to be stiff) in order to yield accurate results, then for all intents and purposes, your problem is stiff. To determine requirements on step size, rely on one of the many libraries out there written by experts (the MATLAB suite is one example, also SUNDIALS, VODE, DASPK, DASSL, LSODE, etc.), which have adaptive time stepping heuristics. The SUNDIALS manual explains the decision rules they use to determine the size of the time steps that package takes, to give you an example of rules that are used in practice.

2. Again, I'd use a library with adaptive time stepping in practice, because it's more efficient to do so. However, if you were coding up a method by yourself, using fixed step sizes, if you noticed large oscillations, or your solution "blowing up", then you'd suspect that your time step was too large, and reduce it. Repeat until you get a reasonably well-behaved numerical solution. Textbooks like Ascher and Petzold and Hairer and Wanner have good examples of this phenomenon.

Answer 2

A better way to look at it is that for a stiff problem, any stable explicit calculation leads to an error that is much smaller than the required error tolerance.

There are many good methods for automatically detecting stiffness using explicit schemes, especially embedded Runge-Kutta pairs. See for example:

• Detecting stiffness with the Fehlberg (4, 5) formulas
• Stiffness detection and estimation of dominant spectrum with explicit Runge-Kutta methods
• Stiffness Detection Strategy for Explicit Runge Kutta Methods

In faleichik's second example, as the step size is reduced, one would see a sudden dramatic decrease in the error to levels far below a typical desired tolerance as the stable timestep threshold is crossed. So a good error estimator would indeed reveal the stiffness of the problem. In the first problem, the error obtained with a stable step size would be in the range of the typical desired tolerance, indicating non-stiffness.

Note as a consequence that any problem becomes non-stiff if a sufficiently strict error tolerance is required.

Answer 3

1. Can we numerically detect stiffness just by applying explicit methods?

• Suppose you have an initial value problem for some ODE on $[0,10]$. You take considerably large stepsize $\tau=1$ and an explicit Euler method, make your calculations with constant step size $\tau$ and get these points:

  

  You estimate the error and it appears to be big. Okay, then you take $\tau=0.1$ and obtain

  The error estimate is acceptable now. Stepsize $\tau=0.1$ is small relatively to $[0,10]^\star$.

  So, is the problem stiff? The answer is NO! Small stepsize here is required to correctly reproduce the oscillations of solution.

  The problem we've solved is $$ y'(t)=-2\cos \pi t,\quad y(0)=1. $$

---

• Now you take another ODE on the same interval, $\tau=1$ and explicit Euler gives you almost the same numerical solution:

  

  You take $\tau=0.1$ and now the numerical solution is

  

  The error estimate now is small. Stepsize $\tau=0.1$ is small relatively to $[0,10]^\star$.

  Is this problem stiff? YES! We've made very much small steps to reproduce the solution which is changing very slowly. This is irrational! The magnitude of time step here is limited by the stability properties of explicit Euler.

  This problem is

  $$ y'(t)=-2 y(t)+\sin t/2,\quad y(0)=1. $$

---

$^\star$ Note that the number of step sizes can be much greater if we take longer interval of integration.

---

Conclusion: the information about timesteps and corresponding errors is not sufficient to detect stiffness. You should also look at the obtained solution. If it varies slowly and stepsize is very small, the problem is most likely to be stiff. If the solution oscillates rapidly and you trust your error estimation technique then this problem is not stiff.

---

2. How to determine the maximum stepsize which allows to integrate stiff problem with explicit method?

If you use some black-box explicit solver with automatic step control then you need nothing to do: the software will take the required stepsize adaptively.

But suppose you want to get the solution manually with constant stepsize, or just want to estimate how many hours you should wait until your explicit method crunches the problem. Then you should know the spectrum of your Jacobi matrix. Suppose it is real and lies in $[\Lambda,0]$ (in your example $\Lambda=-1000$).

Then you should calculate the real stability interval (stability domain in complex case) of your explicit method. It is not too difficult, you need just consult any textbook on this topic. In the case of explicit Euler this interval is $[-2,0]$. Now, if you want your solution not to blow-up, you should take $\tau$ such that $\Lambda\tau$ lies in the stability interval, i. e. in our case $$ \tau\leq \frac{2}{|\Lambda|}. $$

If you want more consistency, you should take $$ \tau\leq \frac{1}{|\Lambda|}, $$ since for $1/|\Lambda|<\tau \leq 2/|\Lambda|$ your solution will most likely produce unnatural (but fading) oscillations.

Of course such an analysis is mostly applicable for linear problems with known spectrum. For more practical problems we should rely on numerical methods of stiffness detection (see references and comments in other answers).