In computational mathematics, a stiff equation is an initial value problem
where , requiring dedicated implicit time stepping methods for its efficient numerical integration. The simplest mathematical characterization of a stiff equation is the necessary condition
Since , where is the Jacobian matrix of at the point , the criterion above is easily evaluated and quantifies stiffness. The criterion is derived, explained and illustrated for nonlinear stiff equations below.
For a linear system with constant coefficients , the divergence is constant, making stiffness a global characteristic whose magnitude is related to the time scale .
For a nonlinear system, stiffness usually varies in space and time along the solution trajectory , where the criterion quantifies stiffness locally. In practical computations, stiff equations are invariably solved using adaptive methods. G Söderlind, LO Jay, M Calvo (2015). "Stiffness 1952--2012: Sixty years in search of a definition." BIT 55, pp 531-558.G Söderlind (2024). "Logarithmic Norms." Berlin-Heidelberg-New York: Springer Series in Computational Mathematics SCM vol 63.
Background
There is a rich literature on stiff differential equations, but intuitive descriptions and heuristics are far more common than attempts at a rigorous definition of the concept. Hairer and WannerE. Hairer, G. Wanner (1996). "Solving orindary differential equations II. Stiff and differential-algebraic problems", 2nd ed. Berlin-Heidelberg-New York: Springer concisely describe the most obvious feature:
Stiff equations are problems for which explicit methods don't work.
This refers to the observation that explicit integration methods are forced to use exceedingly small time steps to maintain numerical stability, preventing such methods from being competitive. Although each step is inexpensive, the total number of steps becomes prohibitively large, and the integration over effectively stalls.
By contrast, implicit methods for stiff equations require costly "algebraic" equation solving on each step. The extra work per step is offset by superior stability properties, allowing the use of large time steps. Without crippling stability restrictions, the total computational effort is manageable, and the requested accuracy can be achieved without loss of efficiency. For some stiff equations, the efficiency may be several orders of magnitude higher than that of even the best explicit methods.
Analogous efficiency issues in other fields of scientific computing
The phenomenon of stiffness is analogous to well-known performance issues in other areas of numerical analysis. For example, in optimization, gradient methods have similar limitations. Using the steepest descent method to minimize a convex functional leads to the iteration
where is the line search step size. This is merely the explicit Euler method for the differential equation . Problems with steep gradients are known for their slow convergence, due to restrictions on the choice of . The typical remedy is to replace the gradient method by some Newton-type method, which essentially corresponds to using an "implicit" method for differential equations.
Another example is in solving nonlinear equations of the form . Simple fixed point iteration, , does not converge for problems where the Lipschitz constant is large, but only for contractions. Yet the equation may have a unique solution under much weaker but common conditions, e.g. when the map is monotone and the logarithmic Lipschitz constant of satisfies . Since fixed point iterations diverge, it is necessary to use a Newton-type iteration.
A third example is iterative solvers for the (discretized) Poisson equation . Here, the Jacobi method (which attempts to avoid matrix algebra) is equivalent to using the explicit Euler method (in pseudo time) for the corresponding method-of-lines discretization of the diffusion equation , seeking its stationary solution. The diffusion equation is a prototypical example of a stiff equation: the time step of the explicit recursion is severely restricted by a CFL condition and required to be proportional to the square of the mesh width in space. This makes the convergence of the Jacobi iteration painfully slow. As before the remedy is to use some form of algebraic equation solving, e.g. by using a preconditioner or a multigrid method.
Thus there are many different mathematical problems, where computational efficiency necessitates the use of dedicated numerical methods, typically of a considerable complexity. Stiff equations is the particular instance for initial value problems in ordinary differential equations. Nevertheless, with a well designed adaptive stiff solver, the efficient numerical integration of most stiff equations is now a routine task.
History and origins of a mathematical characterization
The first mention of stiff equations is found in Curtiss and Hirschfelder in 1952 C.F. Curtiss, J.O. Hirschfelder (1952). "Integration of stiff equations". Proc. Nat. Acad. Sci. 38, pp 235-243. . The authors discuss a scalar model equation equivalent to
and characterize stiffness as follows:
If is the desired resolution of or the interval which will be used in the numerical integration, the equation is "stiff" if
and is well behaved.
This characterization relates a time scale to the decay rate of transients, governed by the coefficient , which is assumed to be negative. As varies in time and space, it allows stiffness to vary accordingly. Prothero and Robinson introduced a similar but more instructive model problem,A. Prothero, A. Robinson (1974). "On the stability and accuracy of one-step methods for solving stiff systems of ordinary differential equations". Math. Comp. 28, 145-162. essentially taking constant and studying the equation
The solution consists of a decaying transient , together with the particular solution , which is supposed to be bounded and smooth. The mathematical solution is then obviously bounded. Contrasting the explicit and implicit Euler methods, one can investigate under what circumstances the numerical methods replicate the properties of the mathematical solution. The Curtiss-Hirschfelder criterion will be considered for and , corresponding to the conditions and , respectively.
The explicit Euler method for a differential equation is defined by
while the implicit Euler method is defined by the recursion
Here is the (constant) time step, and approximates the exact solution at time . In the explicit method is computed directly by an evaluation of the vector field . In the implicit method, however, one must solve the "algebraic" equation for . This improves the stability of the recursion.
For the Prothero-Robinson problem the explicit Euler method produces the numerical solution (on integrated form)
The impediment is clearly visible: positive powers of the factor produce exponential growth (numerical instability) unless we impose the stability condition . Since is real and negative, this implies . Thus the explicit Euler method cannot use step sizes adapted to the regularity of . Instead the total number of steps to complete the integration, , becomes extremely large. Efficiency is inversely proportional to stiffness as quantified by the product .
As a consequence the explicit Euler method may take "forever" to solve a scalar problem with a perfectly smooth mathematical solution once the (fast) transient has decayed. This is also the case even if the initial value is . Thus, the step size restriction is not a matter of whether the transient is present. It is due to the fact that any step size such that inevitably excites numerical instability.
If instead the implicit Euler method is used, the numerical solution becomes
Here, the stability requirement is , which is satisfied for all when , no matter how large is. A minor issue is that in practice one needs to resolve the transient (this may initially call for a short step size), but once it has decayed, one can use "large" step sizes, in the sense that is adapted to the regularity of , irrespective of the magnitude of . Of the two methods above, only the implicit Euler method can handle stiff equations efficiently.
Example Consider the Prothero-Robinson test problem
with exact solution . Taking , the homogeneous solution is not present in the exact solution , but the exponential term will show up in the local solutions passing through the points generated by the numerical methods. Taking the problem is solved on , using and steps, respectively ( and ). The experiment demonstrates the impact of varying , which affects the stability/damping of the numerical methods. We take , using the same computational setup for the explicit and implicit Euler methods.
The results are shown in Figure 1. With , there is a moderately strong exponential damping, as is evident from the fast contracting local solutions. The setup is deliberately understated for visualization purposes; in real stiff equations much larger values would be typical. Nevertheless, is sufficiently large and negative to demonstrate the step size limitations of the explicit Euler method. For , its computed solution becomes oscillatory and diverges from the exact solution, even though the initial value was taken on the exact solution. Repeating the same exercise for , say, shows how the difference between the two methods grows with increasing stiffness.
In more realistic problems, the instability is often more dramatic and referred to as the solution "blowing up". Here, the maximum stable step size () has only been exceeded by a small factor. Meanwhile, the implicit Euler method proceeds without loss of stability, and is able to produce accurate results also for the larger step size.
In real computations, stiff equations are solved using adaptive methods. These select the step size automatically to keep the computational process in check, usually controlling both accuracy and stability. Thus, if an adaptive explicit method is used, the solution may not go unstable, as the adaptivity takes measures not to exceed the maximum stable step size. Instead, the method proceeds at the (prohibitive) cost of low efficiency.
The essence of the analysis is that a problem parameter, , which is characteristic of the vector field and related to the transient behavior of solutions, is linked to a time scale, either the step size or the range of integration directly. The products or are both dimensionless, scaling invariant, and negative. A general stiffness concept must reflect these properties, and also be able to address stiffness in nonlinear differential equations.
Systems of stiff equations
The scalar Prothero-Robinson problem can be extended to linear systems of equations,
as well as nonlinear systems,
provided that . The general stiffness concept needs to specify what is meant by a vector field being "large and negative". The analysis above also shows that the smooth function at most plays a minor role. Thus, by substituting for , we have three types of problems,
Here stiffness needs to be characterized for matrices as well as for nonlinear maps . The first equation above is the well known Dahlquist linear test equation, primarily used to determine the stability region of a time stepping method. The stability region of a method is the set of all such that the numerical method produces bounded solutions. The previous discussion of stiffness is effectively based on the Dahlquist test equation.