Stiff equation


In mathematics, a stiff equation is a differential equation for which certain numerical methods for solving the equation are numerically unstable, unless the step size is taken to be extremely small. It has proven difficult to formulate a precise definition of stiffness, but the main idea is that the equation includes some terms that can lead to rapid variation in the solution.
When integrating a differential equation numerically, one would expect the requisite step size to be relatively small in a region where the solution curve displays much variation and to be relatively large where the solution curve straightens out to approach a line with slope nearly zero. For some problems this is not the case. In order for a numerical method to give a reliable solution to the differential system sometimes the step size is required to be at an unacceptably small level in a region where the solution curve is very smooth. The phenomenon is known as stiffness. In some cases there may be two different problems with the same solution, yet one is not stiff and the other is. The phenomenon cannot therefore be a property of the exact solution, since this is the same for both problems, and must be a property of the differential system itself. Such systems are thus known as stiff systems.

Motivating example

Consider the initial value problem
The exact solution is
We seek a numerical solution that exhibits the same behavior.
The figure illustrates the numerical issues for various numerical integrators applied on the equation.
One of the most prominent examples of the stiff Ordinary differential equations is a system that describes the chemical reaction of Robertson:
If one treats this system on a short interval, for example, there is no problem in numerical integration. However, if the interval is very large, then many standard codes fail to integrate it correctly.
Additional examples are the sets of ODEs resulting from the temporal integration of large chemical reaction mechanisms. Here, the stiffness arises from the coexistence of very slow and very fast reactions. To solve them, the software packages KPP and Autochem can be used.

Stiffness ratio

Consider the linear constant coefficient inhomogeneous system
where and is a constant, diagonalizable, matrix with eigenvalues and corresponding eigenvectors. The general solution of takes the form
where the κt are arbitrary constants and is a particular integral. Now let us suppose that
which implies that each of the terms
as, so that the solution
approaches asymptotically as ;
the term will decay monotonically if λt is real and sinusoidally if λt is complex.
Interpreting x to be time, is called the transient solution and the steady-state solution.
If is large, then the corresponding
term will decay quickly as
x increases and is thus called a fast transient; if
is small, the corresponding term
decays slowly and is
called a slow transient. Let be defined by
so that is the fastest
transient and the
slowest. We now define the stiffness ratio as

Characterization of stiffness

In this section we consider various aspects of the phenomenon of stiffness. "Phenomenon" is probably a more appropriate word than "property", since the latter rather implies that stiffness can be defined in precise mathematical terms; it turns out not to be possible to do this in a satisfactory manner, even for the restricted class of linear constant coefficient systems. We shall also see several qualitative statements that can be made in an attempt to encapsulate the notion of stiffness, and state what is probably the most satisfactory of these as a "definition" of stiffness.
J. D. Lambert defines stiffness as follows:
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.

There are other characteristics which are exhibited by many examples of stiff problems, but for each there are counterexamples, so these characteristics do not make good definitions of stiffness. Nonetheless, definitions based upon these characteristics are in common use by some authors and are good clues as to the presence of stiffness. Lambert refers to these as "statements" rather than definitions, for the aforementioned reasons. A few of these are:
  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.
  3. Stiffness occurs when some components of the solution decay much more rapidly than others.

    Etymology

The origin of the term "stiffness" has not been clearly established. According to Joseph Oakland Hirschfelder, the term "stiff" is used because such systems correspond to tight coupling between the driver and driven in servomechanisms.
According to Richard. L. Burden and J. Douglas Faires,
Significant difficulties can occur when standard numerical techniques are applied to approximate the solution of a differential equation when the exact solution contains terms of the form eλt, where λ is a complex number with negative real part.
...
Problems involving rapidly decaying transient solutions occur naturally in a wide variety of applications, including the study of spring and damping systems, the analysis of control systems, and problems in chemical kinetics. These are all examples of a class of problems called stiff systems of differential equations, due to their application in analyzing the motion of spring and mass systems having large spring constants.

For example, the initial value problem
with m = 1, c = 1001, k = 1000, can be written in the form with n = 2 and
and has eigenvalues
. Both eigenvalues have negative real part and the stiffness ratio is
which is fairly large. System then certainly satisfies statements 1 and 3. Here the spring constant k is large and the damping constant c is even larger.
The exact solution to is
Note that behaves quite nearly as a simple exponential x0e−t, but the presence of the e−1000t term, even with a small coefficient is enough to make the numerical computation very sensitive to step size. Stable integration of requires a very small step size until well into the smooth part of the solution curve, resulting in an error much smaller than required for accuracy. Thus the system also satisfies statement 2 and Lambert's definition.

A-stability

The behaviour of numerical methods on stiff problems can be analyzed by applying these methods to the test equation subject to the initial condition with. The solution of this equation is. This solution approaches zero as when If the numerical method also exhibits this behaviour, then the method is said to be A-stable. A-stable methods do not exhibit the instability problems as described in the motivating example.

Runge–Kutta methods

s applied to the test equation take the form, and, by induction,. The function is called the stability function. Thus, the condition that as is equivalent to. This motivates the definition of the region of absolute stability, which is the set. The method is A-stable if the region of absolute stability contains the set, that is, the left half plane.

Example: The Euler methods

Consider the Euler methods above. The explicit Euler method applied to the test equation is
Hence, with. The region of absolute stability for this method is thus which is the disk depicted on the right. The Euler method is not A-stable.
The motivating example had. The value of z when taking step size is, which is outside the stability region. Indeed, the numerical results do not converge to zero. However, with step size, we have which is just inside the stability region and the numerical results converge to zero, albeit rather slowly.

Example: Trapezoidal method

Consider the trapezoidal method
when applied to the test equation, is
Solving for yields
Thus, the stability function is
and the region of absolute stability is
This region contains the left-half plane, so the trapezoidal method is A-stable. In fact, the stability region is identical to the left-half plane, and thus the numerical solution of converges to zero if and only if the exact solution does. Nevertheless, the trapezoidal method does not have perfect behavior: it does damp all decaying components, but rapidly decaying components are damped only very mildly, because as. This led to the concept of L-stability: a method is L-stable if it is A-stable and as. The trapezoidal method is A-stable but not L-stable. The implicit Euler method is an example of an L-stable method.

General theory

The stability function of a Runge–Kutta method with coefficients and is given by
where denotes the vector with ones. This is a rational function.
Explicit Runge–Kutta methods have a strictly lower triangular coefficient matrix and thus, their stability function is a polynomial. It follows that explicit Runge–Kutta methods cannot be A-stable.
The stability function of implicit Runge–Kutta methods is often analyzed using order stars. The order star for a method with stability function is defined to be the set. A method is A-stable if and only if its stability function has no poles in the left-hand plane and its order star contains no purely imaginary numbers.

Multistep methods

s have the form
Applied to the test equation, they become
which can be simplified to
where z = hk. This is a linear recurrence relation. The method is A-stable if all solutions of the recurrence relation converge to zero when Re z < 0. The characteristic polynomial is
All solutions converge to zero for a given value of z if all solutions w of Φ = 0 lie in the unit circle.
The region of absolute stability for a multistep method of the above form is then the set of all for which all w such that Φ = 0 satisfy |w| < 1. Again, if this set contains the left-half plane, the multi-step method is said to be A-stable.

Example: The second-order Adams–Bashforth method

Let us determine the region of absolute stability for the two-step Adams–Bashforth method
The characteristic polynomial is
which has roots
thus the region of absolute stability is
This region is shown on the right. It does not include all the left half-plane so the Adams–Bashforth method is not A-stable.

General theory

Explicit multistep methods can never be A-stable, just like explicit Runge–Kutta methods. Implicit multistep methods can only be A-stable if their order is at most 2. The latter result is known as the second Dahlquist barrier; it restricts the usefulness of linear multistep methods for stiff equations. An example of a second-order A-stable method is the trapezoidal rule mentioned above, which can also be considered as a linear multistep method.