MarketStiff equation
Company Profile

Stiff equation

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

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 are forced to use exceedingly small time steps h to maintain numerical stability, preventing such methods from being competitive. Although each step is inexpensive, the total number of steps N = T/h becomes prohibitively large, and the integration over [0,T] 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 ==
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 F(x) leads to the iteration x^{k+1} = x^{k} - h\,{\mathrm {grad}}_xF(x^{k}) \, , where h is the line search step size. This is merely the explicit Euler method for the differential equation \dot x = -{\mathrm {grad}}_xF(x^{k}). Problems with steep gradients are known for their slow convergence, due to restrictions on the choice of h. 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 x = g(x). Simple fixed point iteration, x^{k+1} = g(x^k), does not converge for problems where the Lipschitz constant L[g] 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 g-I is monotone and the logarithmic Lipschitz constant of g satisfies M[g] . 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 -\Delta u = f. 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 u_t = \Delta u + f, 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. == Motivating example ==
Motivating example
Consider the initial value problem The exact solution (shown in cyan) is {{NumBlk||y(t)=e^{-15t}, \quad y(t) \to 0 \text{ as } t\to\infty .|}} 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. {{ordered list {{NumBlk||y_{n+1} = y_n+\tfrac12 h\bigl(f(t_n,y_n)+f(t_{n+1},y_{n+1})\bigr), |}} where y'=f(t,y). 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. }} == Another example ==
Another example
One of the most prominent examples of the stiff ordinary differential equations (ODEs) is a system that describes the chemical reaction of Robertson: {{NumBlk||\begin{align} \mathrm{A} \xrightarrow{0.04}& \mathrm{B} \\ \mathrm{B} + \mathrm{B} \xrightarrow{3\cdot 10^7}& \mathrm{C} + \mathrm{B} \\ \mathrm{B} + \mathrm{C} \xrightarrow{1\cdot 10^4}& \mathrm{A} + \mathrm{C} \end{align} \begin{align} \dot x &= - 0.04 x + 10^4 y \cdot z \\ \dot y &= \hphantom{-} 0.04 x - 10^4 y \cdot z - 3\cdot 10^7 y^2 \\ \dot z &= 3\cdot 10^7 y^2 \end{align} If one treats this system on a short interval, for example, t \in [0, 40] 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 ==
Stiffness ratio
Consider the linear constant coefficient inhomogeneous system where \mathbf y, \mathbf f \in \mathbb{R}^n and \mathbf A is a constant, diagonalizable, n \times n matrix with eigenvalues \lambda_t \in \mathbb{C}, t = 1, 2, \ldots , n (assumed distinct) and corresponding eigenvectors \mathbf c_t \in \mathbb{C}^n, t = 1, 2, \ldots , n . The general solution of () takes the form {{NumBlk|| \mathbf y(x) = \sum_{t=1}^{n} \kappa_t e^{\lambda_t x} \mathbf c_t + \mathbf g(x), |}} where the \kappa_t are arbitrary constants and \mathbf g(x) is a particular integral. Now let us suppose that {{NumBlk|| \operatorname{Re}(\lambda_t) |}} which implies that each of the terms e^{\lambda_t x} \mathbf c_t \to 0 as x \to \infty , so that the solution \mathbf y(x) approaches \mathbf g(x) asymptotically as x \to \infty ; the term e^{\lambda_t x} \mathbf c_t will decay monotonically if \lambda_t is real and sinusoidally if \lambda_t is complex. Interpreting x to be time (as it often is in physical problems), \sum_{t=1}^n \kappa_t e^{\lambda_t x} \mathbf c_t is called the transient solution and \mathbf g(x) the steady-state solution. If \left|\operatorname{Re}(\lambda_t)\right| is large, then the corresponding term \kappa_t e^{\lambda_t x} \mathbf c_t will decay quickly as x increases and is thus called a fast transient; if \left| \operatorname{Re}(\lambda_t) \right| is small, the corresponding term \kappa_t e^{\lambda_t x} \mathbf c_t decays slowly and is called a slow transient. Let \overline{\lambda}, \underline{\lambda} \in \{ \lambda_t, t = 1, 2, \ldots , n \} be defined by {{NumBlk|| \left| \operatorname{Re}( \overline{\lambda} ) \right| \geq \left| \operatorname{Re}( \lambda_t ) \right| \geq \left| \operatorname{Re}( \underline{\lambda} ) \right|, \qquad t = 1, 2, \ldots , n |}} so that \kappa_t e^{\overline{\lambda} x} \mathbf c_t is the fastest transient and \kappa_t e^{\underline{\lambda} x} \mathbf c_t the slowest. We now define the stiffness ratio as {{NumBlk|| \frac{ \left| \operatorname{Re}( \overline{\lambda} ) \right| } { \left| \operatorname{Re}( \underline{\lambda} ) \right| }. |}} == Characterization of stiffness ==
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: • A linear constant coefficient system is stiff if all of its eigenvalues have negative real part and the stiffness ratio is large. • Stiffness occurs when stability requirements, rather than those of accuracy, constrain the step length. • Stiffness occurs when some components of the solution decay much more rapidly than others. == Etymology ==
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^{\lambda t}, where \lambda 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). For example, the initial value problem with m = 1, c = 1001, k = 1000, can be written in the form () with n = 2 and {{NumBlk|| \mathbf A = \begin{pmatrix} 0 & 1 \\ -1000 & -1001 \end{pmatrix} , \qquad \mathbf f(t) = \begin{pmatrix} 0 \\ 0 \end{pmatrix} , \qquad \mathbf x(0) = \begin{pmatrix} x_0 \\ 0 \end{pmatrix} , |}} and has eigenvalues \overline{\lambda} = -1000, \underline{\lambda} = -1 . Both eigenvalues have negative real part and the stiffness ratio is {{NumBlk|| \frac{ \left| -1000 \right| }{ \left| -1 \right| } = 1000, |}} 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. (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 () is {{NumBlk|| x(t) = x_0 \left( - \frac{1}{999} e^{-1000 t} + \frac{1000}{999} e^{-t} \right) \approx x_0 e^{-t}. |}} Equation behaves quite similarly to a simple exponential x_0 e^{-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 ==
A-stability
The behaviour of numerical methods on stiff problems can be analyzed by applying these methods to the test equation y' = ky subject to the initial condition y(0) = 1 with k \in \mathbb{C}. The solution of this equation is y(t) = e^{kt}. This solution approaches zero as t\to\infty when \operatorname{Re}(k) If the numerical method also exhibits this behaviour (for a fixed step size), then the method is said to be A-stable. 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
Runge–Kutta methods applied to the test equation y' = k\cdot y take the form y_{n+1} = \phi(hk)\cdot y_n, and, by induction, y_n = \bigl(\phi(hk)\bigr)^n\cdot y_0. The function \phi is called the stability function. Thus, the condition that y_n \to 0 as n \to \infty is equivalent to |\phi(hk)| . This motivates the definition of the region of absolute stability (sometimes referred to simply as stability region), which is the set \bigl\{z \in \mathbb{C} \,\big|\,|\phi(z)|. The method is A-stable if the region of absolute stability contains the set \bigl\{ z \in \Complex \,\big|\, \operatorname{Re}(z) , that is, the left half plane. Example: The Euler methods Consider the Euler methods above. The explicit Euler method applied to the test equation y' = k\cdot y is y_{n+1} = y_n + h\cdot f(t_n, y_n) = y_n + h\cdot(ky_n) = y_n + h\cdot k\cdot y_n = (1+h\cdot k) y_n. Hence, y_n = (1 + hk)^n\cdot y_0 with \phi(z) = 1 + z. The region of absolute stability for this method is thus \bigl\{ z \in \mathbb{C} \,\big|\, |1+z| 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 = \tfrac14 is z = -15\times\tfrac14 = -3.75, which is outside the stability region. Indeed, the numerical results do not converge to zero. However, with step size h = \tfrac18, we have z = -1.875 which is just inside the stability region and the numerical results converge to zero, albeit rather slowly. Example: Trapezoidal method Consider the trapezoidal method y_{n+1} = y_n + \tfrac{1}{2}h\cdot \bigl(f(t_n,y_n)+f(t_{n+1},y_{n+1})\bigr), when applied to the test equation y' = k\cdot y, is y_{n+1} = y_n + \tfrac{1}{2}h\cdot \left(ky_n+ky_{n+1}\right). Solving for y_{n+1} yields y_{n+1}=\frac{1+\frac{1}{2}hk}{1-\frac{1}{2}hk}\cdot y_n. Thus, the stability function is \phi(z)=\frac{1+\frac12 z}{1-\frac12 z} and the region of absolute stability is \left\{ z \in \mathbb{C} \ \left|\ \left| \frac{1+\frac12 z}{1-\frac12 z} \right| 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' = k\cdot y 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 \phi(z) \to 1 as z \to -\infty . This led to the concept of L-stability: a method is L-stable if it is A-stable and |\phi(z)| \to 0 as z \to \infty . 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 \mathbf A and \mathbf b is given by \phi(z) = \frac{\det\left(\mathbf{I} - z\mathbf{A} + z\mathbf{e} \mathbf{b}^\mathsf{T}\right)}{\det \left(\mathbf{I} - z\mathbf{A}\right)}, where \mathbf e 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 \mathbf 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 \bigl\{ z \in \Complex \,\big|\, |\phi(z)| > |e^z| \bigr\}. 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 ==
Multistep methods
Linear multistep methods have the form y_{n+1}=\sum_{i=0}^s a_i y_{n-i} + h \sum_{j=-1}^s b_j f\left(t_{n-j},y_{n-j}\right). Applied to the test equation, they become y_{n+1}=\sum_{i=0}^s a_i y_{n-i} + h k \sum_{j=-1}^s b_jy_{n-j}, which can be simplified to \left(1-b_{-1}z\right)y_{n+1}-\sum_{j=0}^s \left(a_j+b_jz\right)y_{n-j}=0 where 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 \operatorname{Re}(z) . The characteristic polynomial is \Phi(z, w) = w^{s+1}-\sum_{i=0}^s a_iw^{s-i} - z\sum_{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 \in \mathbb{C} for which all w such that \Phi(z,w) = 0 satisfy |w| . 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( \tfrac32 f(t_n, y_n) - \tfrac12 f(t_{n-1}, y_{n-1}) \right) . The characteristic polynomial is \Phi(w,z) = w^2 - \left(1+\tfrac32z\right)w + \tfrac12 z = 0 which has roots w = \tfrac12 \left(1 + \tfrac32 z \pm\sqrt{1 + z + \tfrac94 z^2}\right), thus the region of absolute stability is \left\{ z \in \mathbb{C} \ \left| \ \left| \tfrac12 \left(1 + \tfrac32 z \pm\sqrt{1 + z + \tfrac94 z^2}\right) \right| 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 -1 \le z \le 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 also==
tickerdossier.comtickerdossier.substack.com