Numerical Analysis Lecture (III): Initial Value Problems and Stiffness Part I
It is best to read Numerical Analysis Lecture (II): Numerical Integration first. This Part I covers ordinary differential equation (ODE) initial value problems, basic one-step methods, consistency, stability, and convergence.
3.1 Introduction
Many models in science, engineering, and economics lead to initial value problems for ordinary differential equations.
Initial value problem
Given a function
\[f:[a,b]\times\mathbb{R}^n\to\mathbb{R}^n\]and an initial value $y_0\in\mathbb{R}^n$, find a function
\[y:[a,b]\to\mathbb{R}^n\]whose derivative $y’$ satisfies the ordinary differential equation
\[y'(t)=f(t,y(t)),\qquad t\in[a,b],\]and the initial condition $y(a)=y_0$. We use the English abbreviation IVP (initial value problem) and write
\[\text{(IVP)}\qquad y'(t)=f(t,y(t)),\qquad t\in[a,b], \tag{3.1}\] \[y(a)=y_0. \tag{3.2}\]In many applications, $t$ represents time, so the name “initial value problem” is natural.
Applications
Equations of motion, reaction kinetics, circuit simulation, and many other models.
The following theorem gives a basic existence and uniqueness result for (IVP).
Theorem 3.1.1 (existence and uniqueness)
Let $f:[a,b]\times\mathbb{R}^n\to\mathbb{R}^n$ be continuous. Suppose there is a constant $L>0$ such that
for all $t\in[a,b]$ and all $y,z\in\mathbb{R}^n$ (the Lipschitz condition). Then:
a) By the Picard-Lindelof theorem, for each $y_0\in\mathbb{R}^n$, the (IVP) has a unique solution
\[y\in C^1([a,b];\mathbb{R}^n).\]b) If $y,z$ are solutions with initial values $y(a)=y_0$ and $z(a)=z_0$, then
\[\|y(t)-z(t)\| \le e^{L(t-a)}\|y_0-z_0\|, \qquad \forall t\in[a,b]. \tag{3.3}\]For proofs, see Heuser [3] or Walter [7]. Part b) is a consequence of Gronwall’s lemma.
Remark: Part b) says that the solution depends continuously on the initial value $y_0$.
3.1.1 Basic Concepts of Numerical Methods
To solve (IVP) numerically, divide $[a,b]$ into subintervals:
\[t_j=a+jh,\qquad j=0,1,\ldots,N,\qquad h=\frac{b-a}{N}.\]Integrating (IVP) and writing
\[y_j=y(t_j),\qquad j=0,\ldots,N,\]gives
\[y_{j+1} =y_j+\int_{t_j}^{t_{j+1}}y'(t)\,dt =y_j+\int_{t_j}^{t_{j+1}}f(t,y(t))\,dt. \tag{3.4}\]The integral on the right cannot be evaluated exactly because $y(t)$ is unknown. We therefore approximate it by a quadrature rule and compute approximations
\[u_j\approx y(t_j),\qquad j=1,\ldots,N,\qquad u_0=y_0.\]The error
\[e_j=y(t_j)-u_j\]is called the discretization error.
3.1.2 Some Important Methods
Explicit Euler method
Using the left endpoint rectangle rule in (3.4),
\[\int_{t_j}^{t_{j+1}}f(t,y(t))\,dt \approx h f(t_j,y_j),\]gives the explicit Euler method:
\[u_0:=y_0,\] \[u_{j+1}:=u_j+h f(t_j,u_j), \qquad j=0,\ldots,N-1. \tag{3.5}\]Equivalently, the difference quotient
\[\frac{y(t_{j+1})-y(t_j)}{h}\]is used as an approximation to $y’(t)$ and therefore to $f(t_j,y_j)$.
Explicit Euler is cheap, but it uses only one tangent at the left endpoint. For large step sizes, the error can accumulate quickly.
Implicit Euler method
Using the right endpoint as the rectangle support point gives the implicit Euler method:
\[u_0:=y_0,\] \[u_{j+1}:=u_j+h f(t_{j+1},u_{j+1}), \qquad j=0,\ldots,N-1. \tag{3.6}\]At every step, $u_{j+1}$ must first be determined from this equation.
The slope is evaluated at the right endpoint, which usually improves stability. The price is that one often has to solve a nonlinear equation at each step.
The trapezoidal rule applied to (3.4) gives
\[u_{j+1} =u_j+\frac{h}{2} \left(f(t_j,u_j)+f(t_{j+1},u_{j+1})\right).\]It is implicit because the right-hand side depends on $u_{j+1}$.
The trapezoidal rule averages the slopes at the two endpoints. It is more symmetric than a single Euler step, but it is still implicit.
Heun method, the first second-order Runge-Kutta method
Replacing the unknown right endpoint in the trapezoidal rule by an explicit Euler prediction gives Heun’s method:
\[u_0=y_0,\] \[u_{j+1} =u_j+\frac{h}{2} \left( f(t_j,u_j) +f(t_{j+1},u_j+h f(t_j,u_j)) \right), \qquad j=0,\ldots,N-1.\]Equivalently,
\[u_{j+1}=u_j+\frac{h}{2}(k_1+k_2),\]where
\[k_1=f(t_j,u_j), \qquad k_2=f(t_{j+1},u_j+hk_1).\]The red line starts at the exact point and is drawn as the starting tangent. The blue line is a vector-field sample at the red predicted point; it is not a tangent to the black exact curve.
Improved Euler method (second-order Runge-Kutta method)
Using a midpoint rectangle rule and an Euler half-step prediction gives
\[u_{j+1} =u_j+h f\left(t_j+\frac{h}{2},\, u_j+\frac{h}{2}f(t_j,u_j)\right), \qquad j=0,\ldots,N-1.\]Equivalently,
\[u_{j+1}=u_j+hk_2,\]where
\[k_1=f(t_j,u_j), \qquad k_2=f\left(t_j+\frac{h}{2},u_j+\frac{h}{2}k_1\right).\]The green line represents the vector-field direction sampled at the predicted midpoint. It should not be read as a tangent to the black exact curve.
Classical fourth-order Runge-Kutta method (RK4, fourth-order Runge-Kutta method)
Applying Simpson’s rule and replacing the midpoint and endpoint values by suitable Taylor expansions leads to the widely used classical fourth-order Runge-Kutta method, abbreviated below as RK4:
\[u_{j+1} =u_j+\frac{h}{6}(k_1+2k_2+2k_3+k_4), \qquad j=0,\ldots,N-1,\]where
\[k_1=f(t_j,u_j),\] \[k_2=f\left(t_j+\frac{h}{2},u_j+\frac{h}{2}k_1\right),\] \[k_3=f\left(t_j+\frac{h}{2},u_j+\frac{h}{2}k_2\right),\] \[k_4=f(t_{j+1},u_j+hk_3).\]RK4 gets its accuracy from several vector-field samples. $k_1$ starts from the tangent at the initial point; $k_2,k_3,k_4$ are sampled at predicted points.
Example: Consider a series circuit with inductance $L$, resistance $R$, and capacitance $C$. The current satisfies
\[LC\,I''(t)+RC\,I'(t)+I(t)=0.\]To apply one-step methods, convert this second-order equation into a first-order system by setting
\[y_1(t)=I(t), \qquad y_2(t)=I'(t).\]Then
\[y'(t) = \begin{pmatrix} y_1'(t)\\ y_2'(t) \end{pmatrix} = \begin{pmatrix} y_2(t)\\ -\frac{R}{L}y_2(t)-\frac{1}{LC}y_1(t) \end{pmatrix}.\]Now take
\[L=C=1,\qquad R=\frac14,\]with initial values
\[y_1(0)=I(0)=1,\qquad y_2(0)=I'(0)=1.\]Figure 3.1: Solution and approximation of the oscillator equation
\[I''(t)+\frac14 I'(t)+I(t)=0\]with initial values $I(0)=I’(0)=1$; the left plot uses $n=50$, and the right plot uses $n=100$.
The organization follows Figure 3.1 in the source notes: $n=50$ on the left and $n=100$ on the right. Heun and RK4 follow the exact oscillation much better than explicit Euler.
The exact solution has the form
\[I(t) =e^{-\frac18 t} \left( \cos(\omega t)+\frac{3\sqrt{7}}{7}\sin(\omega t) \right), \qquad \omega=\frac{3\sqrt{7}}{8}.\]Figure 3.1 shows that Heun’s method and RK4 give good approximations, while explicit Euler has a visibly larger error.
3.1.3 Convergence and Consistency
The methods above can be written in the general one-step form
\[u_0=y_0,\] \[u_{j+1} =u_j+h\varphi(t_j,h;u_j,u_{j+1}), \qquad j=0,\ldots,N-1. \tag{3.7}\]Definition 3.1.2
The function $\varphi(t,h;u,v)$ is called the method function in (3.7). If $\varphi$ does not depend on $v$, the method is explicit; otherwise it is implicit.
The German original writes this notation as Die Funktion $\varphi(t,h;u,v)$.
Here $u$ represents the state at the left endpoint of the current step, corresponding to $u_j$ in (3.7); $v$ represents the state at the right endpoint, corresponding to $u_{j+1}$. Explicit methods do not really use $v$, while implicit methods do, so $u_{j+1}$ must be solved from an equation.
The quantity
\[\tau(t,h) = \frac{1}{h} \left( y(t+h)-y(t) -h\varphi(t,h;y(t),y(t+h)) \right)\]is called the local truncation error or consistency error of method (3.7) at position $t$.
Definition 3.1.3
If there are constants $C>0$ and $\bar h>0$ such that
for all $0<h\le\bar h$ and all $t\in[a,b-h]$, then method (3.7) is said to have consistency order $p$ for (IVP).
Consistency asks what happens in one step if the method starts from the exact solution.
The method is called stable if there exists a constant $K>0$ such that
\[\|\varphi(t,h;u,v)-\varphi(t,h;\tilde u,\tilde v)\| \le K\bigl(\|u-\tilde u\|+\|v-\tilde v\|\bigr)\]for all $t\in[a,b]$ and all $u,v,\tilde u,\tilde v\in\mathbb{R}^n$.
Stability here does not mean zero one-step error. It means that small changes in $u$ and $v$ do not cause uncontrolled changes in the method function.
The method has convergence order $p$ if there are constants $M>0$ and $H>0$ such that
\[\|e_j\| =\|y(t_j)-u_j\| \le Mh^p\]for $j=0,\ldots,N$ and all
\[h=\frac{b-a}{N}\le H.\]Convergence concerns the global discretization error at all grid points. If a method converges, reducing $h$ makes the whole numerical trajectory approach the exact one.
Example 3.1.4 (explicit Euler method): Euler’s method has consistency order $1$.
Indeed, if $f\in C^1([a,b]\times\mathbb{R}^n;\mathbb{R}^n)$ and $y$ solves $y’=f(t,y)$, then $y\in C^2([a,b];\mathbb{R}^n)$. Taylor expansion gives
\[y(t+h) =y(t)+f(t,y(t))h +\frac12\bigl(y_i''(t+\xi_i h)\bigr)_{1\le i\le n}h^2.\]Therefore
\[\|\tau(t,h)\| = \left\| \frac1h\bigl(y(t+h)-y(t)-h f(t,y(t))\bigr) \right\| \le Ch.\]3.1.4 A Convergence Theorem
We now state a basic convergence theorem for one-step methods.
Theorem 3.1.5
Let
be the solution of (IVP). Suppose method (3.7) has consistency order $p$ and is stable. Then the method has convergence order $p$. More precisely, there exists $H>0$ such that
\[\|e_j\| =\|y(t_j)-u_j\| \le \frac{e^{4K|t_j-a|}-1}{4K}\,2C h^p\]for $j=0,\ldots,N$ and all
\[h=\frac{b-a}{N}\le H.\]Proof (for interested readers): Let
\[y_j=y(t_j),\qquad e_j=y_j-u_j.\]Consistency gives
\[y_{j+1} =y_j+h\varphi(t_j,h;y_j,y_{j+1})+h\tau(t_j,h).\]Subtracting the numerical method and using stability yields
\[\|e_{j+1}\| \le (1+hK)\|e_j\|+hK\|e_{j+1}\|+hCh^p.\]For sufficiently small $H$, this implies
\[\|e_{j+1}\| \le (1+4hK)\|e_j\|+2Ch^{p+1}.\]The discrete Gronwall lemma below, together with $e_0=0$, gives the stated estimate.
Lemma 3.1.6: For $L>0$, $a_j\ge0$, $h_j>0$, and $b\ge0$, if
\[a_{j+1}\le(1+h_jL)a_j+h_jb, \qquad j=0,1,\ldots,n-1,\]then
\[a_j \le \frac{e^{Lt_j}-1}{L}b +e^{Lt_j}a_0, \qquad t_j:=\sum_{i=0}^{j-1}h_i.\]The induction step is
\[\begin{aligned} a_{j+1} &\le (1+h_jL) \left( \frac{e^{Lt_j}-1}{L}b +e^{Lt_j}a_0 \right) +h_jb\\ &\le \left( \frac{e^{L(t_j+h_j)}-1-h_jL}{L} +h_j \right)b +e^{L(t_j+h_j)}a_0\\ &= \frac{e^{Lt_{j+1}}-1}{L}b +e^{Lt_{j+1}}a_0. \end{aligned}\]3.1.5 Explicit Runge-Kutta Methods
The ideas behind RK4 can be generalized to construct methods with high consistency order.
An $r$-stage explicit Runge-Kutta method has the form
\[u_{j+1} =u_j+h\sum_{i=1}^{r}\beta_i k_i, \qquad k_i=f\left(t_j+\gamma_i h,\, u_j+h\sum_{\ell=1}^{i-1}\alpha_{i\ell}k_\ell \right). \tag{3.7}\]The coefficients are often written in a Butcher tableau:
\[\begin{array}{c|cccc} \gamma_1 & 0 & 0 & \cdots & 0\\ \gamma_2 & \alpha_{21} & 0 & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ \gamma_r & \alpha_{r1} & \alpha_{r2} & \cdots & 0\\ \hline & \beta_1 & \beta_2 & \cdots & \beta_r \end{array}\]Explicit Euler, improved Euler, Heun’s method, and RK4 are all examples of this framework. By choosing the number of stages $r$ and the coefficients appropriately, one can construct explicit Runge-Kutta methods of higher consistency order. The order conditions for $p=1,2,3,4$ are systems of algebraic equations in the coefficients; detailed derivations can be found in Deuflhard and Bornemann [1].
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.
- RK4: fourth-order Runge-Kutta method.
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…