Numerical Analysis Lecture (III): Initial Value Problems and Stiffness Part II
It is best to read Numerical Analysis Lecture (III): Initial Value Problems and Stiffness Part I first. This second part continues with stiffness, stability regions, A-stability, and L-stability for ordinary differential equation (ODE) initial value problems (IVPs).
3.2 Stiff Differential Equations
Stiff systems occur in many applications, for instance in chemical reaction models and in semi-discretizations of partial differential equations. They are still initial value problems, but many numerical methods are forced to use an unacceptably small step size $h$ in order to compute a reliable approximation.
Start from an $n$-dimensional initial value problem. We write IVP$_n$ for an $n$-dimensional initial value problem:
\[\text{(IVP}_n\text{)}\qquad y'(t)=f(t,y(t)),\qquad t\in[a,b],\] \[y(a)=y_0,\]where
\[f:[a,b]\times\mathbb{R}^n\to\mathbb{R}^n, \qquad y_0\in\mathbb{R}^n.\]There is no single universally accepted definition of stiffness. Intuitively, a stiff problem contains two time scales at once: one component changes slowly, while another component decays very quickly.
For the linear case we write LIVP$_n$ for a linear initial value problem:
\[\text{(LIVP}_n\text{)}\qquad y'(t)=Ay(t)+c,\qquad t\in[a,b],\] \[y(a)=y_0,\]where $A\in\mathbb{R}^{n\times n}$ and $c\in\mathbb{R}^n$.
Assume that $A$ is diagonalizable and has eigenvalues $\lambda_i$ with corresponding eigenvectors $v_i$. If $y_P$ is a particular solution, then the general solution has the form
\[y(t)=y_H(t)+y_P(t), \qquad y_H(t)=\sum_{i=1}^{n}C_i e^{\lambda_i t}v_i.\]If
\[\operatorname{Re}(\lambda_i)<0,\qquad i=1,\ldots,n,\]then
\[\lim_{t\to\infty} y_H(t)\to 0,\]so all solutions approach $y_P$. Terms with $\operatorname{Re}(\lambda_i)\ll -1$ decay very quickly, while terms with $\operatorname{Re}(\lambda_i)\not\ll -1$ decay much more slowly. In words: if a system has eigenvalues whose real parts are far below $0$, and also eigenvalues whose real parts are close to $0$ but still negative, then fast and slow scales coexist. Such systems are called stiff systems; see Definition 3.2.2.
In a stiff problem, the fast-decaying component quickly disappears from the visible solution. Its eigenvalue still has a real part far below $0$, however, and this can restrict the stable step size of explicit methods.
Example 3.2.1: Consider
\[y'=Ay, \qquad y(0)=y_0:= \begin{pmatrix} C_1+C_2\\ C_1-C_2 \end{pmatrix},\]where $C_1,C_2\in\mathbb{R}$ and
\[A= \begin{pmatrix} -\frac{1001}{2} & \frac{999}{2}\\ \frac{999}{2} & -\frac{1001}{2} \end{pmatrix}.\]The matrix $A$ has eigenvalues
\[\lambda_1=-1,\qquad \lambda_2=-1000\]with eigenvectors
\[v_1=\begin{pmatrix}1\\1\end{pmatrix}, \qquad v_2=\begin{pmatrix}1\\-1\end{pmatrix}.\]Hence the solution is
\[y(t)= C_1\begin{pmatrix}1\\1\end{pmatrix}e^{-t} +C_2\begin{pmatrix}1\\-1\end{pmatrix}e^{-1000t}.\]The second term is practically irrelevant after a very short time. The first term dominates, and it also tends to $0$ as $t\to\infty$. For a suitable integration method, we would like the numerical values $u_j$ to satisfy
\[\lim_{j\to\infty}u_j=0\]without forcing the step size to be extremely small.
For example, applying the explicit Euler method gives
\[u_0=y_0 =C_1\begin{pmatrix}1\\1\end{pmatrix} +C_2\begin{pmatrix}1\\-1\end{pmatrix}\]and then
\[u_1=(I+hA)u_0 =C_1(1+h\lambda_1)\begin{pmatrix}1\\1\end{pmatrix} +C_2(1+h\lambda_2)\begin{pmatrix}1\\-1\end{pmatrix}.\]Induction gives
\[u_j =C_1(1+h\lambda_1)^j\begin{pmatrix}1\\1\end{pmatrix} +C_2(1+h\lambda_2)^j\begin{pmatrix}1\\-1\end{pmatrix}.\]If $C_2\ne0$, then we must choose
\[|1+h\lambda_2|<1,\]that is,
\[-h\lambda_2=1000h<2,\]to guarantee $\lim_{j\to\infty}u_j=0$. A good method should guarantee this, as far as possible, for all $h>0$.
For a multimode system, the step size must place every $q_i=\lambda_i h$ inside the stability region. The fastest decaying mode often determines the largest explicit step size.
Thus, even though the exact solution itself changes only slowly after the initial transient, Euler’s method still needs a very small step size. Such differential equations are called stiff differential equations. The precise definition is not unique; the following definition is one common version.
Definition 3.2.2: If all eigenvalues of $A$ have non-positive real parts, and $A$ has both eigenvalues satisfying $\operatorname{Re}(\lambda_i)\ll -1$ and eigenvalues whose real parts are close to $0$ but still negative, then the initial value problem (LIVP$_n$) is called stiff.
We now turn to numerical methods for stiff differential equations.
To derive a simple model problem, first consider the case $c=0$ in (LIVP$_n$):
\[y'=Ay, \qquad y(0)=y_0. \tag{3.9}\]Assume that $A$ is diagonalizable. Then there is an $M\in\mathbb{R}^{n\times n}$ such that
\[MAM^{-1}=\operatorname{diag}(\lambda_1,\ldots,\lambda_n),\]where $\lambda_1,\ldots,\lambda_n$ are the eigenvalues of $A$. Let $z=My$. Then
\[z'=MAy=MAM^{-1}z =\operatorname{diag}(\lambda_1,\ldots,\lambda_n)z, \qquad z(0)=My_0.\]Thus the components $z_i$ of $z=My$ satisfy
\[z_i'=\lambda_i z_i, \qquad z_i(0)=(My_0)_i. \tag{3.10}\]For stiff differential equations we also have $\operatorname{Re}(\lambda_i)\le0$, with some eigenvalues whose real parts are far below $0$ and others whose real parts are close to $0$ but still negative.
Observation
If a numerical method behaves well for all scalar differential equations of the form (3.10), then it is usually suitable for the original system as well.
This leads to the model problem
\[y'=\lambda y, \qquad \lambda\in\mathbb{C},\quad \operatorname{Re}(\lambda)<0, \qquad y(0)=1. \tag{3.11}\]The exact solution is
\[y(t)=e^{\lambda t},\]and it satisfies
\[\lim_{t\to\infty}y(t)=0.\]Applying a numerical method with step size $h>0$ produces
\[t_j=jh.\]We want the numerical sequence to satisfy
\[\lim_{j\to\infty}u_j=0. \tag{3.12}\]This leads to the following definition.
Definition 3.2.3 (A-stability, absolute stability)
The term A-stability is also commonly called absolute stability.
If a method applied to the model problem (3.11) produces a sequence ${u_j}_{j\in\mathbb{N}_0}$ satisfying
\[\lim_{j\to\infty}u_j=0\]for every $h>0$, then the method is called A-stable, or absolutely stable.
In general, when a numerical method is applied to the model problem (3.11), one obtains
\[u_{j+1}=R(q)u_j, \qquad q:=\lambda h,\]where $R:D\to\mathbb{C}$ and $0\in D\subseteq\mathbb{C}$.
Read this as follows. Once a one-step method is applied to the test equation $y’=\lambda y$, the whole update often reduces to “multiply the old value by one complex number.” That number is $R(q)$, the amplification factor. It depends on the combined variable $q=\lambda h$, not on $\lambda$ and $h$ separately. In other words, the decay rate of the differential equation and the chosen step size jointly determine whether one numerical step damps the current value, preserves its size, or amplifies it.
The set $D$ is the part of the complex plane where $R$ is defined. Writing $0\in D\subseteq\mathbb{C}$ only says that $R$ is meaningful at least near $q=0$; it need not be defined for every complex $q$. For example, implicit Euler has $R(q)=1/(1-q)$, so $q=1$ is excluded. The implicit trapezoidal rule has $R(q)=(1+q/2)/(1-q/2)$, so $q=2$ is excluded. Stability regions are discussed only where the corresponding $R$ is defined.
With this notation, the stability test is immediate:
\[u_j=R(q)^j u_0.\]| If $ | R(q) | \le1$, one step does not enlarge the current value. If $ | R(q) | <1$, repeated steps damp this mode. If $ | R(q) | >1$, the numerical solution can grow even when the exact solution should decay. The stability region $S$ is therefore the set of all values $q=\lambda h$ for which $ | R(q) | \le1$. |
Example 3.2.4: Applying the explicit Euler method to (3.11) gives
\[u_{j+1}=u_j+h\lambda u_j=(1+q)u_j.\]Thus
\[R(q)=1+q,\]and condition (3.12) becomes
\[|1+q|<1.\]This is exactly the open disk centered at $-1$ with radius $1$ in the complex plane.
The interval $[-2,0]$ on the negative real axis is contained in the stability region. For the explicit Euler method, stability on the negative real axis therefore requires
\[\lambda h\in[-2,0],\]or equivalently
\[h\le \frac{2}{|\lambda|}.\]For a stiff problem this can force $h$ to be very small.
Theorem 3.2.5: Let a numerical method be applied to (3.11), and let $R$ be its stability function. Then the method is A-stable if and only if
\[S:=\{q\in\mathbb{C}: |R(q)|<1\}\]satisfies
\[S\supset\{q\in\mathbb{C}:\operatorname{Re}(q)<0\}.\]Definition 3.2.6 (L-stability)
L-stability is the conventional name used in the numerical ODE literature. Unlike the A in A-stability, the L is usually kept as part of the name rather than expanded into one universally used phrase. The condition to remember is that L-stability is stronger than A-stability: the stability function must also tend to $0$ far along the negative real axis. In other words, for large negative values of $q=\lambda h$, the method should not merely avoid amplification; it should strongly damp the rapidly decaying mode.
A-stability prevents growth throughout the left half-plane. L-stability additionally damps very fast modes close to zero in the numerical solution.
3.2.1 Stability Regions of Some Methods
Explicit Euler method
Applying explicit Euler to the model problem (3.11) gives
\[u_{j+1}=u_j+h\lambda u_j=(1+q)u_j.\]Thus the stability function is
\[R(q)=1+q.\]The stability region is
\[S=\{q\in\mathbb{C}: |1+q|<1\}.\]Therefore explicit Euler is not A-stable, for example at $q=-1+2i$.
Remark 3.2.7: In fact, one can prove that no explicit Runge-Kutta method is A-stable.
Implicit Euler method
Implicit Euler gives
\[u_{j+1}=u_j+h\lambda u_{j+1},\]and therefore
\[u_{j+1}=\frac{1}{1-q}u_j.\]Thus
\[R(q)=\frac{1}{1-q}.\]The stability region is
\[S=\{q\in\mathbb{C}: |1-q|>1\}.\]It contains the entire left half-plane. Also,
\[\lim_{\operatorname{Re}(q)\to -\infty}\frac{1}{1-q}=0.\]Thus implicit Euler is A-stable and even L-stable.
Implicit trapezoidal rule
The implicit trapezoidal rule gives
\[u_{j+1} =u_j+\frac{h}{2}\left(\lambda u_j+\lambda u_{j+1}\right),\]so
\[u_{j+1} = \frac{1+\frac{q}{2}}{1-\frac{q}{2}}u_j.\]The stability function is
\[R(q)=\frac{1+\frac{q}{2}}{1-\frac{q}{2}}.\]The stability region contains the left half-plane, so the method is A-stable. However,
\[\lim_{\operatorname{Re}(q)\to -\infty} \frac{1+\frac{q}{2}}{1-\frac{q}{2}} =-1,\]so it is not L-stable.
Each panel is drawn in the complex plane for $q=\lambda h$: the horizontal axis is $\operatorname{Re}(q)$ and the vertical axis is $\operatorname{Im}(q)$. The implicit Euler stability region is not just the left half-plane; it is the exterior of the disk $|1-q|<1$. The trapezoidal rule has the left half-plane $\operatorname{Re}(q)\le0$ as its stability region.
Implicit Runge-Kutta methods
Implicit Runge-Kutta methods are especially useful for stiff differential equations. If the matrix of Butcher coefficients $\alpha_{ij}$ is not strictly lower triangular, then the method is implicit. Its stage equations are
\[k_i = f\left( t+\gamma_i h,\, u+h\sum_{j=1}^{r}\alpha_{ij}k_j \right), \qquad i=1,\ldots,r.\]The step is then
\[u_{j+1}=u_j+h\sum_{i=1}^{r}\beta_i k_i.\]These equations usually have to be solved as a nonlinear system at each step. In return, one can choose the coefficients so that the method is L-stable and has order $p=2r$.
Explicit Runge-Kutta stages can be computed in order. In an implicit Runge-Kutta method, the $k_i$ usually depend on one another, so the full set of stages is solved first.
Let
\[A=(\alpha_{ij})_{i,j=1}^{r},\qquad \beta=(\beta_1,\ldots,\beta_r)^T, \qquad \mathbf{1}=(1,\ldots,1)^T.\]For the model equation (3.11), the stage equations can be written as
\[k=\lambda u_j\mathbf{1}+qAk,\]or
\[(I-qA)k=\lambda u_j\mathbf{1}.\]Thus
\[k=\lambda u_j(I-qA)^{-1}\mathbf{1}.\]The next value is
\[u_{j+1} = u_j+h\beta^T k = \left(1+q\beta^T(I-qA)^{-1}\mathbf{1}\right)u_j.\]Hence
\[R(q) =1+q\,\beta^T(I-qA)^{-1}\mathbf{1} = \frac{\det(I-qA+q\mathbf{1}\beta^T)} {\det(I-qA)}.\]Therefore $R(q)$ is a rational function.
Return to Numerical Analysis Lecture (III): Initial Value Problems and Stiffness Part I.
References
- [1] P. Deuflhard and F. Bornemann. Numerische Mathematik II. de Gruyter, Berlin, 2002. 3.1.5.
- [2] P. Deuflhard and F. Hohmann. Numerische Mathematik I. de Gruyter, Berlin, 2008. 1.2.3.
- [3] H. Heuser. Gewöhnliche Differentialgleichungen. Teubner, Stuttgart, 1989. 3.1.
- [4] R. Plato. Numerische Mathematik kompakt. Vieweg Verlag, Braunschweig, 2000. 1.2.3, 6.3.2.
- [5] J. Stoer. Numerische Mathematik 1. Springer Verlag, Berlin, 1994. 1.2.3, 4.4.2.
- [6] W. Törnig and P. Spellucci. Numerische Mathematik für Ingenieure und Physiker 2. Springer Verlag, Berlin, 1990. 1.2.3.
- [7] W. Walter. Gewöhnliche Differentialgleichungen. Springer, Berlin, 1986. 3.1.
- [8] J. Werner. Numerische Mathematik 2. Vieweg Verlag, Braunschweig, 1992. 6.1.4.
Abbreviations and Terms
- ODE: ordinary differential equation.
- IVP: initial value problem.
- LIVP: linear initial value problem.
- RK4: fourth-order Runge-Kutta method.
- A-stability: A-stability, also called absolute stability.
- L-stability: L-stability; in addition to A-stability, the stability function tends to $0$ far along the negative real axis.
Source, Copyright, and Usage Notes
This article mainly refers to the numerical analysis lecture notes in TU Darmstadt’s open repository: mathe3-script-2011-SoSe.pdf The upstream repository includes an Unlicense notice. This article is published for personal study, translation, and knowledge organization. The English wording, explanatory additions, and remade figures in this article do not represent the original authors or any official position. The personal organization, English text, explanatory notes, and remade figures in this article may be used for non-commercial study, discussion, and citation with attribution and the original link. Since part of this article is based on translation and organization of TU Darmstadt’s public lecture notes, the original material and any materials it may contain should remain subject to the original authors, repository, and license notices. For commercial use, systematic redistribution, publication, or large-scale adaptation, please verify the licensing status of the original material as well. If there are any translation, formula, terminology, or interpretation errors, or if the rights holder believes the material has been used improperly, please contact me and I will correct or remove it promptly.
AI feedback
AnonymousLoading AI feedback…