# The definition of stiff ODE system

Consider an IVP for ODE system $y'=f(x,y)$, $y(x_0)=y_0$. Most commonly this problem is considered stiff when Jacobi matrix $\frac{\partial f}{\partial y}(x_0,y_0)$ has both eigenvalues with very large negative real part and eigenvalues with very small negative real part (I consider only the stable case).

On the other hand, in the case of just one equation, for example Prothero-Robinson equation $y'=\lambda y + g'+\lambda g$, it is called stiff when $\lambda\ll -1$.

So there are two questions:

1. Why small eigenvalues are included in the definition of stiffness for ODE systems? I believe that the presence of only very large negative real parts is quite enough for system to be stiff, because this makes us use small timesteps for explicit methods.

2. Yes, I know that the most common stiff problems (e. g. arising from parabolic PDEs) do have both large and small eigenvalues. So the second question: is there a good natural example of large stiff system without very small eigenvalues (or alternatively with mild ratio $\lambda_{\max}/\lambda_{\min}$)?

OK, let's modify the question. Consider two two-dimensional linear ODE systems: first with eigenvalues {-1000000,-0.00000001} and second with {-1000000,-999999}. As for me, both of them are stiff. But if we consider stiffness ratio definition, the second system is not. The main question: why stiffness ratio is considered at all?

And the second part of question is still important, lets paraphrase it: I'm looking for a "natural" large ODE system with big negative eigenvalues and mild stiffness ratio (not greater than, say, 100).

• Welcome to scicomp.se. Your questions are answered thoroughly on wikipedia: en.m.wikipedia.org/wiki/Stiff_equation – David Ketcheson Jan 20 '12 at 8:57
• I think between the comment by @DavidKetcheson and the several sources I quoted, you'll see that stiffness ratio is just a guideline. It's not perfect; that's why it's not in the definition. It happens to be a characteristic of many, but not all, stiff systems. And as for the second part, I think you'll be hard-pressed to find it unless it has special structure or arises in an application. I gave you an example of such an application where the stiffness ratio is not always large, and I encourage you to look at Hairer and Wanner's book. – Geoff Oxberry Jan 20 '12 at 10:55
• @David: I cannot agree with you. Take for example one-dimensional problem y'=-50(y-cos x). The "eigenvalue" is -50. One can not solve this problem with explicit Euler with stepsizes greater than 2/50. If we replace -50 with -50000 the restriction on the timestep becomes 2/50000. What "units" can we choose here to overcome this barrier? – faleichik Jan 20 '12 at 15:16
• @faleichik The $\cos x$ part of your example fixes the time scale of the "slow manifold" (which is likely the time scale that you are interested in, though it's conceivable that you would be interested in much shorter time scales). I don't believe it is possible to define stiffness without choosing an observational time scale (perhaps implicitly by stating properties that you want to conserve over longer times). The stiffness ratio only quantifies the scale separation between the fastest and slowest time scales of the autonomous system. – Jed Brown Jan 20 '12 at 18:43
• There is a new, better answer to this question in this paper. – David Ketcheson Jan 15 '17 at 8:12

Stiffness involves some separation of scales. In general, if you are interested in phase of the fastest mode in the system, then you have to resolve it and the system is not stiff. But frequently, you are interested in the long-term dynamics of a "slow manifold" rather than the precise rate at which a solution off the slow manifold approaches it.

Chemical reactions and reacting flows are common examples of stiff systems. The van der Pol oscillator is a common benchmark problem for ODE integrators that has a tunable stiffness paramater.

An ocean is another example that is perhaps helpful to visualize. Tsunamis (surface gravity waves) travel at the speed of an airplane and produce complex wave structure, but dissipate over long time scales and are mostly inconsequential to the long-term dynamics of the ocean. Eddies, or the other hand, travel about 100 times slower at quite pedestrian speeds, but cause mixing and transport temperature, salinity, and biogeochemical tracers that are relevant. But the same physics that propagates a surface gravity wave also supports an eddy (a quasi-equilibrium structure), so eddy velocity, path under Coriolis, and rate of dissipation are dependent on the gravity wave speed. This presents an opportunity for a time integration scheme designed for stiff systems to step over the time scale of the gravity wave and only resolve the relevant dynamical time scales. See Mousseau, Knoll, and Reisner (2002) for discussion of this problem with a comparison of splitting and fully implicit time integration schemes.

Note that diffusive processes are usually considered to be stiff because the fastest time scale in the discrete system is mesh-dependent, scaling with $(\Delta x)^2$, but the time scale of the relevant physics is mesh independent. In fact, the fastest time scales for a given mesh represent spatially local relaxation to the slower manifold on which longer spatial scales evolve, so implicit methods can be very accurate even in strong norms despite not resolving the fastest scales.

## Part 1

Small eigenvalues are not included in the definition of stiffness for ODE (initial value problem) systems. There is no satisfying definition of stiffness that I know of, but the best definitions I've come across are:

If a numerical method with a finite region of absolute stability, applied to a system with any initial conditions, is forced to use in a certain interval of integration a step length which is excessively small in relation to the smoothness of the exact solution in that interval, then the system is said to be stiff in that interval. (Lambert, J. D. (1992), Numerical Methods for Ordinary Differential Systems, New York: Wiley.)

An IVP [initial value problem] is stiff in some interval $$[0,b]$$ if the step size needed to maintain stability of the forward Euler method is much smaller than the step size required to represent the solution accurately. (Ascher, U. M. and Petzold, L. P. (1998), Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations, Philadelphia:SIAM.)

Stiff equations are equations where certain implicit methods, in particular BDF, perform better, usually tremendously better, than explicit ones. (C. F. Curtiss & J. O. Hirschfelder (1952): Integration of stiff equations. PNAS, vol. 38, p. 235-243)

The Wikipedia article on stiff equations goes on to attribute the following "statements" to Lambert:

1. A linear constant coefficient system is stiff if all of its eigenvalues have negative real part and the stiffness ratio is large.

2. Stiffness occurs when stability requirements, rather than those of accuracy, constrain the step length. [Note that this "observation" is essentially the definition from Ascher and Petzold.]

3. Stiffness occurs when some components of the solution decay much more rapidly than others.

Each of these observations has counterexamples (though admittedly I couldn't produce one off the top of my head).

## Part 2

Probably the best example I could come up with would be integrating any sort of large combustion reaction system in chemical kinetics under conditions that result in ignition. The system of equations will be stiff until ignition, and then it will no longer be stiff because the system has passed an initial transient. The ratio of largest to smallest eigenvalue should not be large except around the ignition event, though such systems tend to confound stiff integrators unless you set exceedingly strict integration tolerances.

The book by Hairer and Wanner also gives several other examples in its first section (Part IV, section 1) that illustrate many other examples of stiff equations. (Wanner, G., Hairer, E., Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems (2002), Springer.)

Finally, it's worth pointing out the observation of C. W. Gear:

Although it is common to talk about "stiff differential equations," an equation per se is not stiff, a particular initial value problem for that equation may be stiff, in some regions, but the sizes of these regions depend on the initial values and the error tolerance. (C. W. Gear (1982): Automatic detection and treatment of oscillatory and/or stiff ordinary differential equations. In: Numerical integration of differential equations, Lecture notes in Math., Vol. 968, p. 190-206.)

• Dear Geoff, thanks for the tolerance :-) I wanted to keep my question simple, but eventually came upon being considered inexperienced. In fact I know all these definitions, but. – faleichik Jan 20 '12 at 10:22
• 1. Small eigenvalues implicitly act in the definition of stiffness ratio: it is big when the demoninator is small. 2. For one-dimensional linear case the stiffness ratio is always one, even for stiff equations. 3. Do you have any reference for the chemical kinetics problem you've sugested? And 4. I'll try to clarify the question in the comments. – faleichik Jan 20 '12 at 10:28
• You can find chemical mechanisms in CHEMKIN format here. The problems are large enough that input files are necessary, and the equations are set up automatically using a chemistry package. I suggest using the input files in conjunction with the chemistry package Cantera and the ODE/DAE solver suite SUNDIALS, which are both open source. You can then solve such problems in C++ or MATLAB. – Geoff Oxberry Jan 20 '12 at 10:34
• I personally take the Curtiss-Hirschfelder sentence as my working definition of stiffness; if explicit RK or Adams is taking too long to solve your problem, then it's likely stiff. – J. M. Feb 3 '12 at 3:42

In fact Jed Brown has cleared the question for me. What I'm doing now is just putting his words in the context.

1. Both 2d linear ODE systems from above are stiff (i. e. hard to solve with explicit methods) on relative big time intervals (e. g. [0,1]).

2. The linear systems with large stiffness ratio can be considered "more stiff" because most likely one needs to integrate them on large time interval. This is due to slow components corresponding to the smallest eigenvalues: the solution slowly tends to the steady state, and this steady state is usually important to reach.

3. On the other hand, integration of systems with small stiffness ratio on large intervals is not interestng: in this case the steady state is reached very fast and we can just extrapolate it.

Thanks to all for this discussion!

The absolute magnitude of the eigenvalues (in a linear, autonomous problem) alone has no meaning at all; it's an artifact of the units you choose to express the problem in.

The chain of comments is getting out of control, so I'm making this an answer. I'm not going to answer the full question; as I said, see wikipedia or the other answers here. I'm just answering the bit that says

Consider two two-dimensional linear ODE systems: first with eigenvalues {-1000000,-0.00000001} and second with {-1000000,-999999}. As for me, both of them are stiff. But if we consider stiffness ratio definition, the second system is not. The main question: why stiffness ratio is considered at all?

Okay, let's consider an example of the second case:

$$y_1'(t) = -1000000 y_1(t)$$ $$y_2'(t) = -999999 y_2(t)$$

Now let's consider a time variable with different units: $t^* = 1000000\cdot t$. Then calculus reveals that

$$y_1'(t^*) = - y_1(t^*)$$ $$y_2'(t^*) = -0.999999 y_2(t^*)$$

Note 1: I chose a diagonal system to make it totally obvious, but if you try it with another system with those eigenvalues, you'll see the same effect, since multiplying a matrix by a constant multiplies its eigenvalues by the same constant.

Note 2: I'm not even discussing here whether the system is stiff. I'm just pointing out that your proposed definition of stiffness (i.e., any problem with $|\lambda|\gg1$) makes no sense, since it would mean that stiffness depends on the units in which I choose to express the problem.

• David, you have not considered the interval of integration. Let it be [0,1] in the first case. Assuming explicit Euler's stability limitations, maximum allowed step is 2/1000000. So we need to make at least 500 000 steps. When you scale the time, the maximum stepsize increases to 2, but the whole interval of integration becomes 1 000 000 and we hit the minimum of 500 000 steps again. – faleichik Jan 22 '12 at 11:41
• @faleichik Yes, now you've got it. Stiffness has not to do with the absolute magnitude of the eigenvalues but with their size relative to your timescale of interest, as Jed noted above. – David Ketcheson Jan 22 '12 at 13:34