Stiff equation

Last updated

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.

Contents

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

Explicit numerical methods exhibiting instability when integrating a stiff ordinary differential equation StiffEquationNumericalSolvers.svg
Explicit numerical methods exhibiting instability when integrating a stiff ordinary differential equation

Consider the initial value problem

The exact solution (shown in cyan) is

We seek a numerical solution that exhibits the same behavior.

The figure (right) illustrates the numerical issues for various numerical integrators applied on the equation.

  1. Euler's method with a step size of oscillates wildly and quickly exits the range of the graph (shown in red).
  2. Euler's method with half the step size, , produces a solution within the graph boundaries, but oscillates about zero (shown in green).
  3. The trapezoidal method (that is, the two-stage Adams–Moulton method) is given by
    (3)
    where . Applying this method instead of Euler's method gives a much better result (blue). The numerical results decrease monotonically to zero, just as the exact solution does.

One of the most prominent examples of the stiff ordinary differential equations (ODEs) is a system that describes the chemical reaction of Robertson: [1]

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 (1011 say), then many standard codes fail to integrate it correctly.

Stiffness ratio

Consider the linear constant coefficient inhomogeneous system

where and is a constant, diagonalizable, matrix with eigenvalues (assumed distinct) and corresponding eigenvectors . The general solution of ( 5 ) takes the form

where the 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 is real and sinusoidally if is complex.

Interpreting to be time (as it often is in physical problems), is called the transient solution and the steady-state solution. If is large, then the corresponding term will decay quickly as 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 [2]

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 (and mostly have been) 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. [3]

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. [4] 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 , 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 (mathematical stiffness) systems of differential equations, due to their application in analyzing the motion of spring and mass systems having large spring constants (physical stiffness). [5]

For example, the initial value problem

with , , , can be written in the form ( 5 ) with and

and has eigenvalues . Both eigenvalues have negative real part and the stiffness ratio is

which is fairly large. System ( 10 ) then certainly satisfies statements 1 and 3. Here the spring constant is large and the damping constant is even larger. [6] (while "large" is not a clearly-defined term, but the larger the above quantities are, the more pronounced will be the effect of stiffness.) The exact solution to ( 10 ) is

Equation 13 behaves quite similarly to a simple exponential , but the presence of the term, even with a small coefficient, is enough to make the numerical computation very sensitive to step size. Stable integration of ( 10 ) 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 (for a fixed step size), then the method is said to be A-stable. [7] A numerical method that is L-stable (see below) has the stronger property that the solution approaches zero in a single step as the step size goes to infinity. A-stable methods do not exhibit the instability problems as described in the motivating example.

Runge–Kutta methods

Runge–Kutta methods 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 (sometimes referred to simply as stability region), 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

The pink disk shows the stability region for the Euler method. Stability region for Euler method.svg
The pink disk shows the stability region for the Euler method.

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

The pink region is the stability region for the trapezoidal method. Stability region for trapezoidal method.svg
The pink region is the stability region for the 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. [8]

General theory

The stability function of a Runge–Kutta method with coefficients and is given by

where denotes the vector with all ones. This is a rational function (one polynomial divided by another).

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. [9]

Multistep methods

Linear multistep methods have the form

Applied to the test equation, they become

which can be simplified to

where . This is a linear recurrence relation. The method is A-stable if all solutions of the recurrence relation converge to zero when . The characteristic polynomial is

All solutions converge to zero for a given value of if all solutions of 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 such that satisfy . 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

The pink region is the stability region for the second-order Adams-Bashforth method. Stability region for AB2 method.svg
The pink region is the stability region for 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 (in fact it only includes the real axis between ) 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. [10]

See also

Notes

  1. Robertson, H. H. (1966). "The solution of a set of reaction rate equations". Numerical analysis: an introduction. Academic Press. pp. 178–182.
  2. Lambert (1992 , pp. 216–217)
  3. Lambert (1992 , pp. 217–220)
  4. Hirshfelder (1963)
  5. Burden & Faires (1993 , p. 314)
  6. Kreyszig (1972 , pp. 62–68)
  7. This definition is due to Dahlquist (1963).
  8. The definition of L-stability is due to Ehle (1969).
  9. The definition is due to Wanner, Hairer & Nørsett (1978); see also Iserles & Nørsett (1991).
  10. See Dahlquist (1963).

Related Research Articles

<span class="mw-page-title-main">Laplace's equation</span> Second-order partial differential equation

In mathematics and physics, Laplace's equation is a second-order partial differential equation named after Pierre-Simon Laplace, who first studied its properties. This is often written as or where is the Laplace operator, is the divergence operator, is the gradient operator, and is a twice-differentiable real-valued function. The Laplace operator therefore maps a scalar function to another scalar function.

<span class="mw-page-title-main">Navier–Stokes equations</span> Equations describing the motion of viscous fluid substances

The Navier–Stokes equations are partial differential equations which describe the motion of viscous fluid substances. They were named after French engineer and physicist Claude-Louis Navier and the Irish physicist and mathematician George Gabriel Stokes. They were developed over several decades of progressively building the theories, from 1822 (Navier) to 1842–1850 (Stokes).

<span class="mw-page-title-main">Ellipsoid</span> Quadric surface that looks like a deformed sphere

An ellipsoid is a surface that can be obtained from a sphere by deforming it by means of directional scalings, or more generally, of an affine transformation.

In mathematical optimization, the method of Lagrange multipliers is a strategy for finding the local maxima and minima of a function subject to equation constraints. It is named after the mathematician Joseph-Louis Lagrange.

<span class="mw-page-title-main">Heat equation</span> Partial differential equation describing the evolution of temperature in a region

In mathematics and physics, the heat equation is a certain partial differential equation. Solutions of the heat equation are sometimes known as caloric functions. The theory of the heat equation was first developed by Joseph Fourier in 1822 for the purpose of modeling how a quantity such as heat diffuses through a given region.

Linear elasticity is a mathematical model of how solid objects deform and become internally stressed due to prescribed loading conditions. It is a simplification of the more general nonlinear theory of elasticity and a branch of continuum mechanics.

<span class="mw-page-title-main">Green's function</span> Impulse response of an inhomogeneous linear differential operator

In mathematics, a Green's function is the impulse response of an inhomogeneous linear differential operator defined on a domain with specified initial conditions or boundary conditions.

In mathematics, the Schwarzian derivative is an operator similar to the derivative which is invariant under Möbius transformations. Thus, it occurs in the theory of the complex projective line, and in particular, in the theory of modular forms and hypergeometric functions. It plays an important role in the theory of univalent functions, conformal mapping and Teichmüller spaces. It is named after the German mathematician Hermann Schwarz.

Geometrical optics, or ray optics, is a model of optics that describes light propagation in terms of rays. The ray in geometrical optics is an abstraction useful for approximating the paths along which light propagates under certain circumstances.

In mathematics, an Euler–Cauchy equation, or Cauchy–Euler equation, or simply Euler's equation is a linear homogeneous ordinary differential equation with variable coefficients. It is sometimes referred to as an equidimensional equation. Because of its particularly simple equidimensional structure, the differential equation can be solved explicitly.

In numerical analysis, the Crank–Nicolson method is a finite difference method used for numerically solving the heat equation and similar partial differential equations. It is a second-order method in time. It is implicit in time, can be written as an implicit Runge–Kutta method, and it is numerically stable. The method was developed by John Crank and Phyllis Nicolson in the mid 20th century.

Linear multistep methods are used for the numerical solution of ordinary differential equations. Conceptually, a numerical method starts from an initial point and then takes a short step forward in time to find the next solution point. The process continues with subsequent steps to map out the solution. Single-step methods refer to only one previous point and its derivative to determine the current value. Methods such as Runge–Kutta take some intermediate steps to obtain a higher order method, but then discard all previous information before taking a second step. Multistep methods attempt to gain efficiency by keeping and using the information from previous steps rather than discarding it. Consequently, multistep methods refer to several previous points and derivative values. In the case of linear multistep methods, a linear combination of the previous points and derivative values is used.

Numerical continuation is a method of computing approximate solutions of a system of parameterized nonlinear equations,

In mathematics, the spectral theory of ordinary differential equations is the part of spectral theory concerned with the determination of the spectrum and eigenfunction expansion associated with a linear ordinary differential equation. In his dissertation, Hermann Weyl generalized the classical Sturm–Liouville theory on a finite closed interval to second order differential operators with singularities at the endpoints of the interval, possibly semi-infinite or infinite. Unlike the classical case, the spectrum may no longer consist of just a countable set of eigenvalues, but may also contain a continuous part. In this case the eigenfunction expansion involves an integral over the continuous part with respect to a spectral measure, given by the Titchmarsh–Kodaira formula. The theory was put in its final simplified form for singular differential equations of even degree by Kodaira and others, using von Neumann's spectral theorem. It has had important applications in quantum mechanics, operator theory and harmonic analysis on semisimple Lie groups.

A differential equation is a mathematical equation for an unknown function of one or several variables that relates the values of the function itself and its derivatives of various orders. A matrix differential equation contains more than one function stacked into vector form with a matrix relating the functions to their derivatives.

In the finite element method for the numerical solution of elliptic partial differential equations, the stiffness matrix is a matrix that represents the system of linear equations that must be solved in order to ascertain an approximate solution to the differential equation.

In mathematics, the method of steepest descent or saddle-point method is an extension of Laplace's method for approximating an integral, where one deforms a contour integral in the complex plane to pass near a stationary point, in roughly the direction of steepest descent or stationary phase. The saddle-point approximation is used with integrals in the complex plane, whereas Laplace’s method is used with real integrals.

In mathematics, the Neumann–Poincaré operator or Poincaré–Neumann operator, named after Carl Neumann and Henri Poincaré, is a non-self-adjoint compact operator introduced by Poincaré to solve boundary value problems for the Laplacian on bounded domains in Euclidean space. Within the language of potential theory it reduces the partial differential equation to an integral equation on the boundary to which the theory of Fredholm operators can be applied. The theory is particularly simple in two dimensions—the case treated in detail in this article—where it is related to complex function theory, the conjugate Beurling transform or complex Hilbert transform and the Fredholm eigenvalues of bounded planar domains.

In mathematics, a linear recurrence with constant coefficients sets equal to 0 a polynomial that is linear in the various iterates of a variable—that is, in the values of the elements of a sequence. The polynomial's linearity means that each of its terms has degree 0 or 1. A linear recurrence denotes the evolution of some variable over time, with the current time period or discrete moment in time denoted as t, one period earlier denoted as t − 1, one period later as t + 1, etc.

Tau functions are an important ingredient in the modern mathematical theory of integrable systems, and have numerous applications in a variety of other domains. They were originally introduced by Ryogo Hirota in his direct method approach to soliton equations, based on expressing them in an equivalent bilinear form.

References