↑ arvau

Wikipedia / Stiff equation

In computational mathematics, a stiff equation is an initial value problem

, requiring dedicated implicit time stepping methods for its efficient numerical integration. The simplest mathematical characterization of a stiff equation is the necessary condition

is the Jacobian matrix, 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.

For a nonlinear system, where the criterion quantifies stiffness locally. In practical computations, stiff equations are invariably solved using adaptive methods.

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 Wanner 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 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 leads to the iteration

is the line search step size. This is merely the explicit Euler method for the differential equation. 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. Simple fixed point iteration, does not converge for problems where the Lipschitz constant is monotone and the logarithmic Lipschitz constant. 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 . The authors discuss a scalar model equation equivalent to

and characterize stiffness as follows:

or the interval which will be used in the numerical integration, the equation is "stiff" if

is well behaved.

varies in time and space, it allows stiffness to vary accordingly. Prothero and Robinson introduced a similar but more instructive model problem, constant and studying the equation

, respectively.

The explicit Euler method is defined by

while the implicit Euler method is defined by the recursion

is computed directly by an evaluation of the vector field. This improves the stability of the recursion.

For the Prothero-Robinson problem the explicit Euler method produces the numerical solution (on integrated form)

.

. Thus, the step size restriction is not a matter of whether the transient is present inevitably excites numerical instability.

If instead the implicit Euler method is used, the numerical solution becomes

. Of the two methods above, only the implicit Euler method can handle stiff equations efficiently.

Example Consider the Prothero-Robinson test problem

, using the same computational setup for the explicit and implicit Euler methods.

The results are shown in Figure 1., say, shows how the difference between the two methods grows with increasing stiffness.

) 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.

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,

, we have three types of problems,

. The first equation above is the well known Dahlquist linear test equation such that the numerical method produces bounded solutions. The previous discussion of stiffness is effectively based on the Dahlquist test equation.

. The step size has to be selected small enough to fulfill this requirement, see Figure 2.

Yet, it is common to find a different characterization of stiffness, in terms of the "stiffness ratio", defined by

. Although stiff systems often do have a large stiffness ratio, Byrne and Hindmarsh emphasize that this is neither necessary nor sufficient to be a stiff equation.

, which is common in some applications. (An example is chemical reaction kinetics, where the singularity corresponds to conservation of mass.)

Thus the stiffness ratio is misleading: it fails to characterize the most important aspects of stiff equations.

has a large stiffness ratio. Naturally, this too fails, as noted by Artemiev and Averina:

] satisfy the above conditions. However, the famous van der Pol equation, which is frequently used as an example of stiff ODEs, does not satisfy this definition ... it is impossible to determine stiffness only by means of eigenvalues of the Jacobi matrix for nonlinear systems

Although in part valid, this argument too is off the mark, as the zero solution of the van der Pol equation is an unstable equilibrium. Thus stiffness does not occur near the origin, but only along the limit cycle as will be seen below, where the van der Pol equation is analyzed in full detail. A precise stiffness characterization must therefore be local, capturing the behavior of (small) perturbations near the actual solution.

Modern characterization: The stiffness indicator

, called the stiffness indicator. This characterization is based on additional observations. Thus Dekker and Verwer write (original emphasis):

The essence of stiffness is that the solution to be computed is slowly varying but that perturbations exist which are rapidly damped.

Shampine expresses a similar view, but is conceptually more specific, relating stiffness to dissipation, and what in practice is irreversibility:

A way we prefer for describing the latter condition is that [the solution] is very unstable in the opposite direction [of time].

The reverse instability can easily be characterized using logarithmic norms, the norm of the solution can be bounded using differential inequalities. Then

, respectively. These are the extrema of a symmetric quadratic form,

. Its stationary points are found from the symmetric eigenvalue problem

are referred to as the logarithmic values, the logarithmic values satisfy

together with the trace identity

it follows from the differential inequality above that

, while the maximum reverse time.

Söderlind et al. now quantify Shampine's observation. Thus, in a stiff system,

. The stiffness indicator is defined

and a stiff equation is characterized by the condition

is positive, offsetting stiffness. This is then accounted for by defining the stiffness indicator so as to have odd parity.

is extremely small. In other words, the flow is nearly a semigroup is "close to singular". Unlike in a Lipschitz system, where one can solve (theoretically) the differential equation also in reverse time, this is effectively impossible in a stiff system. The reverse time stiff equation corresponds to an extremely ill conditioned "inverse problem".

Using logarithmic norms for nonlinear maps is defined globally, the stiffness indicator is instead defined locally along the trajectory.

In addition to the stiffness indicator, Söderlind et al. introduce a local reference time scale defined by

, referred to as the (local) stiffness factor, but there may also be parts of the integration where the step size stiffness factor is only moderate.

An important aspect of implicit time stepping methods is the need to solve equations of the form

, which represents stiff computation. Rewriting the equation as

. Since

, it may require a minor step size restriction. In the Shampine criterion, as well as in the stiffness indicator, it holds that

is not a contraction), existence and uniqueness are usually guaranteed also when the step size stiffness factor is large.

A simplified stiffness indicator

, an inexpensive alternative is of interest. Noting that the stiffness indicator is the arithmetic average of the largest and the smallest logarithmic value, an option is to replace this average by the average of all logarithmic values. With the aforementioned trace identity (see above), we define the scaled (arithmetic average) real trace functional

is both inexpensive and trivial: no eigenvalue computations are required. Moreover, for a nonlinear map,

serves as an alternative stiffness indicator. Stiffness can therefore also be characterized by the divergence condition

which is readily evaluated for most problems.

, it holds that

. The differential equation for the determinant of the flow can now be rewritten

Hence

, phase volume is compressed.

is a necessary condition for stability being the simplest.

Examples of nonlinear stiff equations and their properties

A large number of nontrivial test problems have been collected in the Bari Test Set for Initial Value Problem Solvers. This contains detailed examples of problems, solvers, and results achieved under specified conditions by well-known codes.

Here three stiff nonlinear problems are selected to illustrate how the (simplified) stiffness indicator is used to identify and characterize stiffness correctly. This includes the ability to distinguish the complex and rapidly changing behavior at (locally unstable) turning points, such as those observed in the van der Pol and Oregonator equations. Further, the stiffness indicator also distinguishes between subintervals where stiffness is prominent and where the problem is nonstiff. Finally in some cases it is even possible to identify which of the dependent variables contributes most to stiffness.

Operational criteria, such as how the adaptive time step depends on the accuracy requirements and varies along the solution, are also illustrated while quantifying the step size stiffness factor accordingly. Taken together, a detailed insight can be derived from the stiffness indicator, some of it even before the actual computation begins.

Problem 1. (Flame propagation) Consider the scalar flame propagation model,

. This simplified model is commonly used as a test problem to challenge the time step adaptivity of software for stiff equations, checking whether the sharp transition between a nonstiff and a stiff region is correctly resolved.

. Since

.

, stability is regained, and the solution enters the stiff regime, see Figure 3 the solution moves.

.

Problem 2. (van der Pol equation) The van der Pol equation is usually presented in a couple of alternate forms, with different scalings. It is a nonlinear problem whose solutions approach a limit cycle. Here we scale time to write the system as

, accounting for a full period. The problem originates in the study of nonlinear oscillations in electric circuits.

(see Figure 4, given by

. The stiffness indicator is readily identified in the second equation of the van der Pol system, hiding in plain sight.

, the latter variable has no particular influence on stiffness.

it follows that stiffness only occurs for, i.e., along the two stable branches, where stiffness is proportional to for a stiff solver.

Note that all this information is available a priori. It is confirmed by the numerical solution of the problem.

Problem 3. (Oregonator equation). The equations are

.

do not coincide, although the differences are small. More importantly, because the equation is a more intricate example of stiffness, we need to gain preliminary insight by keeping the parameters as long as possible in the analysis. We therefore choose the scaled divergence, which offers this advantage. Thus we find

where

will contribute depends on what magnitude it will reach. In Figure 5.

.

When an adaptive third order A-stable Runge-Kutta method is used, the step sizes exceed the reference time scale by three to four orders of magnitude. This corresponds to the efficiency gain of the stiff solver over the potential use of an explicit method.

Methods and software for stiff equations

.

The stability region is determined by applying the method to the Dahlquist test equation contains the entire left half-plane, is called A-stable.

. Therefore one usually has to settle for less. Among multistep methods the best choice for stiff equations are the backward differentiation are regularly used.

Software based on BDF-type methods (or similar) are Matlab's ode15s, and C or Fortran codes such as CVODE, LSODE, MEBDF, DASSL. The software is elaborate, robust and of high complexity, uses variable order and variable time step adaptivity, and usually offers many options for problems with special structure, or for special handling of Jacobian matrices. Accuracy criteria can be specified both for absolute and relative error estimates. Some codes include options for nonstiff equations, or for differential-algebraic equations.

is referred to as the stability function. If

the method is A-stable. There are A-stable methods of high order, but few implementations. An advantage of A-stable Runge-Kutta methods is that

is negative definite, then the method produces a contractive recursion. (This is a variant of the von Neumann inequality in a nonlinear problem, then the method produces a contractive recursion. Other special stability requirements are also common, such as L-stability..

, and the Fortran code RADAU5 order Radau IIa method. For a general survey of diagonally implicit Runge-Kutta methods used for stiff equations, see Kennedy and Carpenter.

In all cases, the practical numerical solution of stiff equations needs dedicated, professional software. Even though it is a routine task, success may often require a careful study of the suitable settings of the codes and a proper understanding of the problem to be solved.

Notes and remarks

1. Is stiffness fully understood? In the literature it is sometimes suggested that no precise definition of stiffness exists. This appears to be due to the widespread mention or use of the stiffness ratio, whose obvious shortcomings are long since recognized. Sometimes operational criteria such as the actual discretization method and accuracy requirements are brought in to resolve the issues. But this over-complicates what is basically a simple issue. Thus, to the practitioner it is obvious that there is a distinction between nonstiff and stiff equations; it must therefore be possible to describe this distinction mathematically. The stiffness indicator provides a simple necessary criterion characterizing and quantifying this distinction. Today it is safe to say that stiffness is a complex but fully understood phenomenon, with excellent, efficient, and reliable dedicated software available.

2. Stiffness indicator.. While the former only involves the two extreme logarithmic values, the latter makes use of the intermediate logarithmic values as well, in spite of the fact that they do not have a similar impact on stability. The advantage of the scaled divergence is that it bypasses eigenvalue computations and supports an a priori analytical understanding of the dynamics, as illustrated in the examples above. An approximate quantification is usually sufficient, as stiffness typically is a matter of orders of magnitude.

3. Dissipative vs. conservative systems. The notion of stiffness as discussed above is a measure of dissipativity, i.e., of damping or energy loss. Some sources suggest that there are also stiff conservative systems, in particular if derived from method-of-lines discretizations of hyperbolic PDEs. But such systems belong in a different category, of (potentially) highly oscillatory systems. With eigenvalues located far out on the imaginary axis, the corresponding high frequencies can be suppressed by choosing an implicit method using large steps. However, in accordance with the sampling theorem, high frequency phenomena are then no longer resolved and wave forms may be distorted. The distinction between dissipative and conservative problems is significant, and the computational approach is different. While stiffness is closely related to the irreversibility of parabolic problems, hyperbolic problems have no intrinsic damping and feature invariants, such as conserving energy. Similarly, separable Hamiltonian systems have a divergence free vector field and hence are "area preserving" (conserve phase volume). Due to this structure, the stiffness indicator is zero.

4. Etymology. The term "stiff" was introduced by Curtiss and Hirschfelder. According to Hirschfelder, the term was chosen because the first examples were associated with servo systems, in which there was a "tight coupling" between the servo and the controlled system. A similar effect can also be seen in control systems with a high gain negative feedback, where higher gain increases stiffness. Another suggested connection is the (mechanical) notion of stiffness, as represented by the second order equation

, the equation is rewritten as

are

is large. For critical damping leads to

Interestingly, computing the stiffness indicator (here identical to the scaled divergence) of the vector field one obtains

which is large if and only if is large, independent of the spring constant. Thus the single cause of stiffness is the damper, which dissipates energy. It is difficult to avoid the conclusion that the term "stiff" is a misnomer, which quickly got established in a new context with a different meaning. While in mechanical engineering a large spring constant typically goes with a large damping constant in order to have (near) critical damping, the confusion can be understood.

See also

Notes

References