To control the global error \( y(t_n)-y_n \) is notoriously difficult, and far beyond what will be discussed in this course. To control the local error in each step and adjust the step size accordingly is rather straightforward, as we will see.
As we already have seen, the local error is found by finding the power series in \( h \) of the exact and the numerical solution. The local error is of order \( p \) if the lowest order terms in the series where the exact and the numerical solution differ is of order \( p+1 \). So the local errors of the two methods are $$ \begin{align*} \mb{y}(t_{n+1};t_n,\mb{y}_n) - \mb{y}_{n+1} &= \mb{\Psi}(t_n,\mb{y}_n)h^{p+1} +\dotsc, \\ \mb{y}(t_{n+1};t_n,\mb{y}_n) - \widehat{\mb{y}}_{n+1} &= \hphantom{\mb{\Psi}(t_n,\mb{y}_n)h^{p+1} } + \dotsc, \end{align*} $$ where \( \mb{\Psi}(t_n,\mb{y}_n) \) is a term consisting of method parameters and differentials of \( \mb{f} \) and \( \dotsc \) contains all the terms of the series of order \( p+2 \) or higher. Taking the difference gives $$ \widehat{\mb{y}}_{n+1} - \mb{y}_{n+1} = \mb{\Psi}(t_n,\mb{y}_n)h^{p+1} + \ldots. $$
Assume now that \( h \) is small, such that the principal error term \( \mb{\Psi}(t_n,\mb{y}_n)h^{p+1} \) dominates the error series. Then a reasonable approximation to the unknown local error \( \mb{l}_{n+1} \) is the local error estimate \( \mb{le}_{n+1} \): $$ \mb{le}_{n+1} = \widehat{\mb{y}}_{n+1} - \mb{y}_{n+1} \approx \mb{y}(t_{n+1};t_n,\mb{y}_n) - \mb{y}_{n+1}. $$
Example: Apply Euler's method of order 1 and Heun's method of order 2 with \( h=0.1 \) to the equation $$ y' = -2ty, \qquad y(0)=1. $$ Use this to find an approximation to the error after one step.
Euler's method: $$ y_1 = 1.0 - 0.1\cdot 2 \cdot 0 \cdot 1.0 = 1.0. $$ Heun's method $$ \begin{align*} k_1 &= -2\cdot 0.0 \cdot 1.0 = 0.0, \\ k_2 &= -2\cdot 0.1\cdot (1+0.0) = -0.2, \\ \widehat{y}_1& = 1.0 + \frac{0.1}{2}\cdot(0.0 - 0.2) = 0.99. \end{align*} $$ The error estimate and the local error are respectively $$ le_{1} = \widehat{y}_1 - y_1 = -10^{-2}, \qquad l_1 = y(0.1)-y_1 = e^{-0.1^2}-1.0 = -0.995 \cdot 10^{-2}. $$ So in this case the error estimate is a quite decent approximation to the actual local error.
Essentially:
Given \( t_n, \mb{y}_n \) and a step size \( h_n \).
From the discussion above, we have that $$ \| \mb{le}_{n+1} \| \approx D h_{n}^{p+1}. $$ where \( \mb{le}_{n+1} \) is the error estimate we can compute, \( D \) is some unknown quantity, which we assume almost constant from one step to the next. We aim for a step size \( h_{\textrm{new}} \) such that $$ \text{Tol} \approx D h_{\textrm{new}}^{p+1}. $$ From these two approximations we get: $$ \frac{\text{Tol}}{\|\mb{le}_{n+1}\|} \approx \left(\frac{h_{\textrm{new}}}{h_n}\right)^{p+1} \qquad \Rightarrow \qquad h_{\textrm{new}} \approx \left( \frac{\text{Tol}}{\|\mb{le}_{n+1}\|} \right)^{\frac{1}{p+1}} h_{n}. $$
To avoid too many rejected steps, it is wise to be a bit conservative when choosing the new step size, so the following is used in practice: $$ h_{\textrm{new}} = P\cdot \left( \frac{\text{Tol}}{\|\mb{le}_{n+1}\|} \right)^{\frac{1}{p+1}} h_{n}. $$ where the pessimist factor \( P < 1 \) is some constant, normally chosen between 0.5 and 0.95.
heun_euler can be written asode_adaptive is written to make it simple to test other pairs of methods. This is also the reason why the function heun_euler returns the order of the lowest order method.using \( \text{Tol}=10^{-3} \) and \( h_0=0.1 \).
See the function ode_example4() in numode.py.
Numerical exercises:
The error estimate is then given by $$ \mb{le}_{n+1} = h\sum_{i=1}^s (\widehat{b}_i - b_i)\mb{k}_i. $$
Using this notation, the Heun-Euler pair is presented as $$ \begin{array}{c|cc} 0 & & \\ 1 & 1 & \\ \hline & 1 & 0 \\ \hline \displaystyle & \frac{1}{2} & \frac{1}{2} \end{array} $$
Let us start this section by an example:
Numerical example 1s: Given the following system of two ODEs $$ \begin{align*} y_1' &= -2y_1 + y_2 + 2\sin(t), & y_1(0) &= 2, \\ y_2' &= (a-1) y_1 - a y_2 + a\,\big(\cos(t)-\sin(t) \big), & y_2(0) &= 3, \end{align*} $$ where \( a \) is some positive parameter. The exact solution, which is independent of the parameter, is $$ y_1(t) = 2 e^{-t} + \sin(t), \qquad y_2(t) = 2e^{-t} + \cos(t). $$ Solve this problem now with some adaptive ODE solver, for instance the Heun-Euler scheme.
Now try the tolerances \( \text{Tol}=10^{-2}, \, 10^{-4}, \, 10^{-6} \), and perform the experiment with two different values of the parameters, \( a=2 \) and \( a=999 \).
See \verb+ode_example_1s+ in \verb+numode.py+.
For \( a=2 \) the expected behaviour is observed, for higher accuracy, more steps are required. But the example \( a=999 \) requires much more steps, and the step size seems almost independent of the tolerance, at least for \( \text{Tol}=10^{-2}, \, 10^{-4} \).
The example above with \( a=999 \) is a typically example of a stiff ODE. What defines these types of ODEs is that there are (at least) two different time scales at play at the same time: a slow time scale that dominates the time evolution of the solution of the ODE, and a fast time scale at which small perturbations of the solution may occur. In physical systems, this might be due to very strong damping of selected components of the system.
Numerical example 1, cont.: The general solution of the ODE can be shown to be $$ \mb{y}(t) = c_1 \left(\begin{array}{c} 1 \\ 1 \end{array}\right) e^{-t} + c_2 \left(\begin{array}{c} -1 \\ a-1 \end{array}\right) e^{-(a+1)t} + \left(\begin{array}{c} \sin(t) \\ \cos(t) \end{array}\right) $$ for some constants \( c_1 \), \( c_2 \). The terms \( e^{-t} \), \( \sin(t) \), and \( \cos(t) \) evolve at a time scale of order \( 1 \). In contrast, the term \( e^{-(a+1)t} \) reverts back to being essentially equal to zero at a time scale of order \( 1/(a+1) \).
When a stiff ODE is solved by some explicit adaptive method like the Heun-Euler scheme, an unreasonably large number of steps is required, and this number seems independent of the tolerance. The problem is that, for explicit methods, the local error is dominated by what is happening at the fast time scale, and the step length will be adapted to that time scale as well. Even worse, any larger step size will lead to instabilities and exponentially increasing oscillations in the numerical solution.
In the remaining part of this note we will explain why this happens, and how we can overcome the problem. For simplicity, the discussion is restricted to linear problems, but also nonlinear ODEs can be stiff, and often will be.
Exercise s1: Repeat the experiment on the Van der Pol equation $$ \begin{align*} y_1' &= y_2, & y_1(0) &= 2, \\ y_2' &= \mu(1-y_1^2)y_2 - y_1, & y_2(0) &= 0. \end{align*} $$ Use \( \mu=2 \), \( \mu=5 \) and \( \mu=50 \).
The analytic solution of this problem is $$ y(x) = y_0\,e^{\lambda x} = y_0 \, e^{\Re \lambda\, x} \bigl(\cos(\Im \lambda\, x) + i \sin(\Im \lambda \, x)\bigr). $$ Since \( \Re\lambda < 0 \), the solution \( y(x) \) tends to zero as \( x\rightarrow \infty \). We want a similar behaviour for the numerical solution, that is $ |y_{n}| \rightarrow 0$ when \( n\rightarrow \infty \).
One step of some Runge–Kutta method applied to the linear test equation can always be written as $$ y_{n+1} = R(z)y_n, \qquad z=\lambda h. $$ The function \( R(z) \colon \mathbb{C} \to \mathbb{C} \) is called the stability function of the method. \( R(z) \) is either a polynomial or a rational function of \( z \).
Example 2s: The application of Euler's method for the linear test equation results in the sequence $$ y_{n+1} = y_n + h \lambda y_n = (1+h\lambda) y_n = (1+z) y_n \qquad\text{ with } z = h\lambda. $$ The stability function of Euler's method is therefore the function $$ R(z) = 1+z. $$
For a comparison, Heun's method for this test equation is $$ \begin{align*} k_1 &= \lambda y_n,\\ k_2 &= \lambda (y_n + h k_1),\\ y_{n+1} &= y_n + \frac{h}{2}(k_1 + k_2), \end{align*} $$ which can be rewritten as $$ y_{n+1} = y_n + \frac{h}{2}(\lambda y_n + \lambda(y_n+h\lambda y_n) = y_n + h\lambda y_n + \frac{(h\lambda)^2}{2} y_n. $$ As a consequence, the stability function for Heun's method is $$ R(z) = 1 + z + \frac{z^2}{2}. $$ One step with the Trapezoidal rule is $$ y_{n+1} = y_n + h\frac12(\lambda y_n + \lambda y_{n+1}) $$ and the corresponding stability function is $$ R(z) = \frac{1+\frac12 z}{1-\frac12 z} $$
We now return back to the analysis of the behaviour of an arbitrary Runge-Kutta method with stability function \( R \). Taking the absolute value on each side of the expression $$ y_{n+1} = R(z) y_n, $$ we see that there are three possible outcomes: $$ \begin{align*} |R(z)| & < 1 \quad \Rightarrow & |y_{n+1}| & < |y_n| \quad \Rightarrow && y_n \rightarrow 0 &&\text{(stable)} \\ |R(z)| &= 1 \quad \Rightarrow & |y_{n+1}| & = |y_n| \\ |R(z)| &> 1 \quad \Rightarrow & |y_{n+1}| &> |y_n| \quad \Rightarrow && |y_n| \rightarrow \infty && (\text{unstable}) \end{align*} $$
The stability region of a method is defined by $$ \mathcal{S} = \{ z \in \mathbb{C} \; :\; |R(z)| \leq 1 \}. $$ To get a stable numerical solution, we have to choose the step size \( h \) such that \( z=\lambda h\in \mathcal{S} \).
Example 2s, continued: In the case of Euler's method, we have obtained the stability function $$ R(z) = 1+z. $$ As a consequence, the stability region for Euler's method is $$ \mathcal{S} = \{ z \in \mathbb{C}\; :\; \lvert 1+z \rvert \leq 1\}. $$ This is a ball in the complex plane, which is centered at \( -1 \) and has a radius of \( 1 \).
We are given a system of \( m \) differential equation of the form $$ \mb{y}' = A \mb{y} + \mb{g}(x). \tag{*} $$ Such systems have been discussed in Mathematics 3, and the technique for finding the exact solution will shortly be repeated here:
Solve the homogenous system \( \mb{y}' = A \mb{y} \), that is, find the eigenvalues \( \lambda_i \) and the corresponding eigenvectors \( \mb{v}_i \) satisfying $$ A\mb{v}_i = \lambda_i \mb{v}_i, \qquad i=1,2,\dotsc,m. \tag{**} $$
Assume that \( A \) has a full set of linearly independent (complex) eigenvectors \( \mb{v}_i \) with corresponding (complex) eigenvalues \( \lambda_i \). Let \( V=[\mb{v}_1,\dots,\mb{v}_m] \), and \( \Lambda = \text{diag}\{\lambda_1,\dotsc,\lambda_m\} \). Then \( V \) is invertible and $$ AV = V\Lambda \qquad\text{ and therefore }\qquad V^{-1}AV = \Lambda. $$
The ODE (*) can thus be rewritten as $$ V^{-1} \mb y' = V^{-1}A V V^{-1} \mb{y}+ V^{-1}\mb{g}(t). $$ Let \( \mb{z} = V^{-1}\mb{y} \) and \( \mb{q}(t)=V^{-1}\mb{g}(t) \) such that the equation can be decoupled into a set of independent scalar differential equations $$ \mb{z}' = \Lambda \mb{z} + \mb{q}(t) \qquad \text{ or, equivalently } \qquad z_i' = \lambda_i z_i + q_i(t), \quad i=1,\dotsc, m. $$ The solution of such equations has been discussed in Mathematics 1. When these solutions are found, the exact solution is given by $$ \mb{y}(t) = V \mb{z}(t), $$ and possible integration constants are given by the initial values.
As it turns out, the eigenvalues \( \lambda_i \in \mathbb{C} \) are the key to understanding the behaviour of the ODE integrators, and it motivates the study of the stability properties of the very simplified, though complex, linear test equation $$ y' = \lambda y, \qquad \lambda \in \mathbb{C}^{-}. $$
The discussion below is also relevant for nonlinear ODEs \( \mb{y}'(t)=\mb{f}(t,\mb{y}(t)) \), in which case we have to consider the eigenvalues of the Jacobian \( \mb{f}_{\mb{y}} \) of \( \mb{f} \) with respect to \( \mb{y} \).
Example 1: We now return to the introductory example. There, the ODE can be written as $$ \mb{y}' = A \mb{y} + \mb{g}(t), $$ with $$ A = \left(\begin{array}{cc} -2 & 1 \\ a-1 & -a \end{array}\right), \qquad \mb{g}(t) = \left(\begin{array}{c}\sin(t) \\ a(\cos(t)-\sin(t)) \end{array}\right). $$ The eigenvalues of the matrix \( A \) are \( \lambda_1 = -1 \) and \( \lambda_2 = -(a+1) \). The general solution is given by $$ \mb{y}(t) = c_1 \left(\begin{array}{c} 1 \\ 1 \end{array}\right) e^{-t} + c_2 \left(\begin{array}{c} -1 \\ a-1 \end{array}\right) e^{-(a+1)t} + \left(\begin{array}{c} \sin(t) \\ \cos(t) \end{array}\right). $$ In the introductory example, the initial values were chosen such that \( c_1=2 \) and \( c_2=0 \). However, for large values of \( a \), the term \( e^{-(a+1)t} \) will still go to 0 almost immediately, even if \( c_2\not=0 \). It is this term that creates problems for the numerical solution.
Numerical example 2: We have already discussed the stability function and stability region for Euler's method in the example above. We now solve the introductory problem $$ \mb{y}' = \left(\begin{array}{cc} -2 & 1 \\ a-1 & -a \end{array}\right) \mb{y} + \left(\begin{array}{c}\sin(t) \\ a(\cos(t)-\sin(t)) \end{array}\right), \qquad \mb{y}(0) = \left(\begin{array}{c} 2 \\ 3 \end{array}\right), \qquad a>0. $$ by Euler's method. We know that the eigenvalues of the matrix \( A \) are \( \lambda_1 = -1 \) and \( \lambda_2 = -(1+a) \).
For the numerical solution to be stable for both eigenvalues, we have to require that the step length \( h \) satisfies $$ \lvert 1 + \lambda_1 h \lvert \le 1 \qquad\text{ and }\qquad \lvert 1 + \lambda_2 h \lvert \le 1. $$ Since both eigenvalues in this case are real and negative, we see after a short computation that this results in the requirement that $$ h \leq \frac{2}{1+a}. $$
Try \( a=9 \) and \( a=999 \). Choose step sizes a little bit over and under the stability boundary, and you can experience that the result is sharp. If \( h \) is just a tiny bit above, you may have to increase the interval of integration to see the unstable solution.
It is the term corresponding to the eigenvalue \( \lambda_2=-(a+1) \) which makes the solution unstable. And the solution oscillate since \( R(z) < -1 \) for \( h>2/(1+a) \).
Exercise 2:
In an ideal world, we would prefer the stability region to satisfy $$\mathcal S \supset \mathbb{C}^- := \{ z \in \mathbb{C}\;:\; \Re z \le 0\},$$ such that the method is stable for all \( \lambda \in \mathbb{C} \) with \( \Re \lambda \le 0 \) and for all \( h \). Such methods are called \( A \)-stable. For all explicit methods, like Euler's and Heun's, the stability function will be a polynomial, and \( |R(z)|\rightarrow \infty \) as \( \Re z\rightarrow -\infty \). Explicit methods can never be \( A \)-stable, and we therefore have to search among implicit methods. The simplest of those is the implicit, or backward, Euler's method, given by $$ y_{n+1} = y_n + hf(t_{n+1}, y_{n+1}). $$ Applied to the linear test equation \( y'=\lambda y \), this results in the update $$ y_{n+1} = y_n + h\lambda y_{n+1} \qquad \text{ or } \qquad y_{n+1} = \frac{1}{1-h\lambda}y_n. $$ We therefore see that we have the stability function $$ R(z) = \frac{1}{1-z}. $$ The stability region for the implicit Euler function is thus $$ \mathcal{S} = \Bigl\{ z \in \mathbb{C} \;:\; \Bigl\lvert \frac{1}{1-z}\Bigr\rvert \le 1\Bigr\} = \{ z \in \mathbb{C}\;:\; \lvert 1-z \rvert \ge 1\}. $$ This is the whole complex plan apart from an open ball with center \( +1 \) and radius \( 1 \). Thus the method is \( A \)-stable, as every complex number \( z \) with \( \Re z \le 0 \) is contained in \( \mathcal{S} \).
In the implementation below, the right hand side of the ODE is implemented as a
function rhs, returning the matrix \( A \) and the vector \( \mb{g}(x) \) for each
step. The function implicit_euler performs one step with implicit Euler. It has
the same interface as the explicit method, so that the function ode_solve can be used as before.
Numerical example 3: Solve the test equation with $$ A = \left(\begin{array}{cc} -2 & 1 \\ a-1 & -a \end{array}\right), \qquad \mb{g}(t) = \left(\begin{array}{c}\sin(t) \\ a(\cos(t)-\sin(t)) \end{array}\right), $$ by the implicit Euler method. Choose \( a=2 \) and \( a=999 \), and try different stepsizes like \( h=0.1 \) and \( h=0.01 \). Are there any stability issues in this case?
Exercise 2: The trapezoidal rule is an implicit method which for a general ODE \( \mb{y}'(x)=f(x,\mb{y}(x)) \) is given by $$ \mb{y}_{n+1} = \mb{y}_{n} + \frac{h}{2}\bigg( \mb{f}(x_n,\mb{y}_n) + \mb{f}(x_{n+1},\mb{y}_{n+1})\bigg). $$
trapezoidal_ieuler.
The interface is as for the embedded pair heun_euler, so the adaptive solver ode_adaptive can be used as before.
Numerical example 4:
Repeat the experiment from the introduction, using trapezoidal_ieuler.
We observe that there are no longer any step size restriction because of stability. The algorithm behaves as expected.
Comment: Implicit methods can of course also be applied for nonlinear ODEs. Implicit Euler's method will be $$ \mb{y}_{n+1} = \mb{y}_n + h\mb{f}(x,\mb{y}_{n+1}), $$ which is a nonlinear system which has to be solved for each step. Similar for the trapezoidal rule. Usually these equations are solved by Newton's method or some simplification of it.