Stiff equation

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.

Motivating example

Consider the initial value problem:,y'(t)=-15y(t),quad t ge 0, y(0)=1.The exact solution is y(t)=e^{-15t}, and clearly y(t) → 0 as t → ∞. We want the numerical solutions to exhibit the same behaviour.

Figure 1 illustrates the numerical issues for various numerical integrators applied on the equation. First, Euler's method with a step size of h=1/4 oscillates wildly and quickly exits the range of the graph. After halving the step size and re-running the integration with h=1/8, the resulting solution is within the graph boundaries, but it oscillates about zero, and by no means represents the exact solution.

The trapezoidal method, that is, the two-stage Adams–Moulton method, is given by:y_{n+1}=y_n+{1over 2}hleft(f(t_n,y_n)+f(t_{n+1},y_{n+1}) ight).Applying this method instead of Euler's method gives a much better result. As the figure shows, the numerical results decrease monotonically to zero, just like the exact solution.

Characterization of stiffness

Certain types of problems can be characterized as stiff:
* problems of the form:: y' = ky + f(t),,:: where k∈C and |k| is large.
* systems of the form:: mathbf{y}' = Kmathbf{y}+mathbf{f}(t),:: where K is a square matrix having at least one eigenvalue l∈C and |l| is large.
* systems of the form:: mathbf{y}' = mathbf{f}(mathbf{y}, t),:: with the Jacobian of f having at least one eigenvalue m∈C and |m| is large.

A-stability

The behaviour of numerical methods on stiff problems can be analyzed by applying these methods to the test equation y' = ky with "k" &isin; C. The solution of this equation is y(t) = mathrm{e}^{kt}. This solution approaches zero as t oinfty when Re k < 0. If the numerical method also exhibits this behaviour, then the method is said to be A-stable. [This definition is due to harvtxt|Dahlquist|1963.] 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 y' = ky take the form y_{n+1}=phi(hk)y_n, and by induction, y_n = phi(hk)^ny_0 (the function &phi; is called the "stability function"). Thus, the condition that y_n &rarr; 0 as n &rarr; &infin; is equivalent to |phi(hk)|<1. This motivates the definition of the "region of absolute stability" (sometimes referred to simply as "stability region"), which is the set {z &isin; C | |phi(z)|<1 }. The method is A-stable if the region of absolute stability contains the set { z &isin; C | Re(z) < 0 }, that is, the left half plane.

Example: The Euler and trapezoidal methods

Consider both the Euler and trapezoidal methods above. The Euler method applied to the test equation y' = ky is: y_{n+1} = y_n + hf(t_n, y_n) = y_n + h(ky_n) = y_n + hky_n = (1+hk)y_n. , Hence, y_n = (1+hk)^ny_0 with phi(z) = 1+z. The region of absolute stability for this method is thus { z &isin; C | |1+z| < 1 } which is the disk depicted on the right. The Euler method is not A-stable.

The motivating example had k=-15. The value of "z" when taking step size h = 1/4 is z = -3.75, which is outside the stability region. Indeed, the numerical results do not converge to zero. However, with step size h = 1/8, we have z = -1.875 which is just inside the stability region and the numerical results converge to zero, albeit rather slowly.

The trapezoidal method: y_{n+1} = y_n + frac12h left(f(t_n,y_n)+f(t_{n+1},y_{n+1}) ight), when applied to the test equation y' = ky, is: y_{n+1} = y_n + frac12h left(ky_n+ky_{n+1}) ight). Solving for y_{n+1} yields:y_{n+1}={1+{1over 2}hk over 1-{1over 2}hk}y_n. Thus, the stability function is:phi(z)={1+{1over 2}z over 1-{1over 2}z}and the region of absolute stability is:left{ z in Bbb{C} left| left| {1+{1over 2}z over 1-{1over 2}z} ight| < 1 ight. ight}.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 y' = ky converges to zero if "and only if" the exact solution does. Nevertheless, the trapezoidal method does not have perfect stability: it does damp all decaying components, but rapidly decaying components are damped only very mildly, because phi(z) o 1 as z o -infty . This led to the concept of L-stability: a method is L-stable if it is A-stable and phi(z) o 0 as |z| o infty . The trapezoidal method is not L-stable. The implicit Euler method is an example of an L-stable method. [The definition of L-stability is due to harvtxt|Ehle|1969.]

General theory

The stability function of a Runge–Kutta method with coefficients "A" and "b" is given by: phi(z) = frac{det(I-zA+zeb^T)}{det(I-zA)}, where "e" denotes the vector with ones. This is a rational function (one polynomial divided by another).

Explicit Runge–Kutta methods have a strictly lower triangular coefficient matrix "A" 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 phi is defined to be the set {z &isin; C | |phi(z)| > mathrm{e}^z }. 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. [The definition is due to harvtxt|Wanner|Hairer|Nørsett|1978; see also harvtxt|Iserles|Nørsett|1991.]

Multistep methods

Linear multistep methods have the form:y_{n+1}=sum_{i=0}^s a_i y_{n-i}+hsum_{j=-1}^s b_j f(t_{n-j},y_{n-j}).Applied to the test equation, they become:y_{n+1}=sum_{i=0}^s a_i y_{n-i}+hksum_{j=-1}^s b_jy_{n-j},which can be simplified to:(1-b_{-1}z)y_{n+1}-sum_{j=0}^s (a_j+b_jz)y_{n-j}=0where z=hk. This is a linear recurrence relation. The method is A-stable if all solutions {y_n} of the recurrence relation converge to zero when Re "z" < 0. The characteristic polynomial is :Phi(z, w) = w^{s+1}-sum_{i=0}^s a_iw^{s-j} - zsum_{j=-1}^s b_jw^{s-j}.All solutions converge to zero for a given value of z if all solutions w of Phi(z,w) = 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 z &isin; C for which all w such that Phi(z,w)=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: y_{n+1} = y_n + h left( frac32 f(t_n, y_n) - frac12 f(t_{n-1}, y_{n-1}) ight) .The characteristic polynomial is: Phi(w,z) = w^2 - (1+ frac32z)w + frac12 z = 0which has roots: w = frac12 Big(1 + frac32 z pmsqrt{1 + z + frac94 z^2}Big), thus the region of absolute stability is: left{ z in Bbb{C} left| left| frac12 Big(1 + frac32 z pmsqrt{1 + z + frac94 z^2}Big) ight| < 1 ight. ight}.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 z=-1 and z=0) 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. [See harvtxt|Dahlquist|1963.]

ee also

* Explicit and implicit methods

Notes

References

*.
*.
*.
*.
*.


Wikimedia Foundation. 2010.

Игры ⚽ Поможем решить контрольную работу

Look at other dictionaries:

  • Stiff — may refer to:*Stiffness, a material s resistance to bending *Stiff differential equation, an equation that exhibits behaviour at two widely different scales (the differential equations describing stiff materials are stiff differential equations)… …   Wikipedia

  • Ordinary differential equation — In mathematics, an ordinary differential equation (or ODE) is a relation that contains functions of only one independent variable, and one or more of their derivatives with respect to that variable. A simple example is Newton s second law of… …   Wikipedia

  • Differential algebraic equation — In mathematics, differential algebraic equations (DAEs) are a general form of (systems of) differential equations for vector–valued functions x in one independent variable t, where is a vector of dependent variables and the system has as many… …   Wikipedia

  • Numerical ordinary differential equations — Illustration of numerical integration for the differential equation y = y,y(0) = 1. Blue: the Euler method, green: the midpoint method, red: the exact solution, y = et. The step size is h = 1.0 …   Wikipedia

  • List of numerical analysis topics — This is a list of numerical analysis topics, by Wikipedia page. Contents 1 General 2 Error 3 Elementary and special functions 4 Numerical linear algebra …   Wikipedia

  • Euler method — In mathematics and computational science, the Euler method, named after Leonhard Euler, is a first order numerical procedure for solving ordinary differential equations (ODEs) with a given initial value. It is the most basic kind of explicit… …   Wikipedia

  • List of mathematics articles (S) — NOTOC S S duality S matrix S plane S transform S unit S.O.S. Mathematics SA subgroup Saccheri quadrilateral Sacks spiral Sacred geometry Saddle node bifurcation Saddle point Saddle surface Sadleirian Professor of Pure Mathematics Safe prime Safe… …   Wikipedia

  • Euler equations (fluid dynamics) — In fluid dynamics, the Euler equations govern inviscid flow. They correspond to the Navier Stokes equations with zero viscosity and heat conduction terms. They are usually written in the conservation form shown below to emphasize that they… …   Wikipedia

  • Numerical stability — In the mathematical subfield of numerical analysis, numerical stability is a desirable property of numerical algorithms. The precise definition of stability depends on the context, but it is related to the accuracy of the algorithm. A related… …   Wikipedia

  • Quark-Stern — Ein Quarkstern, auch Seltsamer genannt, ist ein theoretisch möglicher Endzustand der Sternentwicklung vor einem Schwarzen Loch. Mit dem Verbrauch seines nuklearen Brennmaterials (Kernfusion) wird die Materie eines Sterns durch die Gravitation… …   Deutsch Wikipedia

Share the article and excerpts

Direct link
Do a right-click on the link above
and select “Copy Link”