Solving ordinary differential equations numerically
Cours ENPC — Pratique du calcul scientifique
Introduction
Objective
$$
% % %
%
$$
\[ \definecolor{gris_C}{RGB}{96,96,96} \definecolor{blanc_C}{RGB}{255,255,255} \definecolor{pistache_C}{RGB}{176,204,78} \definecolor{pistache2_C}{RGB}{150,171,91} \definecolor{jaune_C}{RGB}{253,235,125} \definecolor{jaune2_C}{RGB}{247,208,92} \definecolor{orange_C}{RGB}{244,157,84} \definecolor{orange2_C}{RGB}{239,119,87} \definecolor{bleu_C}{RGB}{126,151,206} \definecolor{bleu2_C}{RGB}{90,113,180} \definecolor{vert_C}{RGB}{96,180,103} \definecolor{vert2_C}{RGB}{68,141,96} \definecolor{pistache_light_C}{RGB}{235,241,212} \definecolor{jaune_light_C}{RGB}{254,250,224} \definecolor{orange_light_C}{RGB}{251,230,214} \definecolor{bleu_light_C}{RGB}{222,229,241} \definecolor{vert_light_C}{RGB}{216,236,218} \]
Goal of this chapter
Study numerical methods to solve initial value problems: \[ \left \{ \begin{aligned} & \mathbf{\boldsymbol{x}}'(t) = \mathbf{\boldsymbol{f}}\bigl(t, \mathbf{\boldsymbol{x}}(t)\bigr), \\ & \mathbf{\boldsymbol{x}}(t_0) = \mathbf{\boldsymbol{x}}_0. \end{aligned} \right. \qquad(1)\] where \(\mathbf{\boldsymbol{f}}\colon \mathbf R\times \mathbf R^n \to \mathbf R^n\).
Application examples
Lagrangian mechanics;
Celestial mechanics;
Electric circuits;
Time-dependent partial differential equations: \[ \left \{ \begin{aligned} \partial_t u &= \partial_x^2 u \\ u(0, x) &= u_0(x). \end{aligned} \right. \xrightarrow{\text{discretization}} \left \{ \begin{aligned} \mathbf{\boldsymbol{x}}' &= \mathsf A \mathbf{\boldsymbol{x}}. \\ \mathbf{\boldsymbol{x}}(0) &= \mathbf{\boldsymbol{x}}_0. \end{aligned} \right. \] \({\color{green} \rightarrow}\) ubiquitous in science.
Remark: higher-order differential equations
Initial value problems of higher order can be recast in form Equation 1.
Example: charged particle (wih mass and charge 1) in electric field: \[ \mathbf{\boldsymbol{x}}''(t) = \mathbf{\boldsymbol{E}}\bigl(\mathbf{\boldsymbol{x}}(t)\bigr). \] Introducing \(\mathbf{\boldsymbol{v}} = \mathbf{\boldsymbol{x}}'\), we rewrite this equation as \[ \left \{ \begin{aligned} & \mathbf{\boldsymbol{v}}'(t) = \mathbf{\boldsymbol{E}}\bigl(\mathbf{\boldsymbol{x}}(t)\bigr), \\ & \mathbf{\boldsymbol{x}}'(t) = \mathbf{\boldsymbol{v}}(t). \end{aligned} \right. \]
Remark: a useful particular case
When \(\mathbf{\boldsymbol{f}}(t, \mathbf{\boldsymbol{x}}) = \mathbf{\boldsymbol{f}}(t)\), the solution is an integral: \[ \mathbf{\boldsymbol{x}}(t) = \mathbf{\boldsymbol{x}}_0 + \int_{t_0}^{t} \mathbf{\boldsymbol{f}}(t) \, \mathrm dt. \] \({\color{green} \rightarrow}\) useful to guess order of convergence of numerical methods
Analysis of the continuous problem
Definition: local, maximal, global solution
A function \(\mathbf{\boldsymbol{x}} \colon I \to \mathbf R^n\) is a local solution if \(t_0 \in I\) and \[ \left \{ \begin{aligned} & \mathbf{\boldsymbol{x}}'(t) = \mathbf{\boldsymbol{f}}\bigl(t, \mathbf{\boldsymbol{x}}(t)\bigr) \qquad \forall t \in I \\ & \mathbf{\boldsymbol{x}}(t_0) = \mathbf{\boldsymbol{x}}_0. \end{aligned} \right. \] A local solution \(\mathbf{\boldsymbol{x}} \colon I \to \mathbf R^n\) is maximal if it has no extension to a larger interval, and global if \(I = \mathbf R\).
Example
Consider the equation \[ \left \{ \begin{aligned} & x'(t) = x(t)^2, \\ & x(0) = 1. \end{aligned} \right. \] The solution \(x_* \colon (-\infty, 1) \to \mathbf R\) given by \[ \mathbf{\boldsymbol{x}}_*(t) = \frac{1}{1-t}. \] is maximal but not global.
Theorem: existence and uniqueness
Suppose that
The function \(\mathbf{\boldsymbol{f}}\) is continuous;
For any bounded set \(D \subset \mathbf R\times \mathbf R^n\) containing \((t_0, \mathbf{\boldsymbol{x}}_0)\), there is \(L(D)\) such that \[ \forall \bigl((t, \mathbf{\boldsymbol{x}}_1), (t, \mathbf{\boldsymbol{x}}_2)\bigr) \in D \times D, \qquad {\lVert { \mathbf{\boldsymbol{f}}(t, \mathbf{\boldsymbol{x}}_1) - \mathbf{\boldsymbol{f}}(t, \mathbf{\boldsymbol{x}}_2) } \rVert} \leqslant L {\lVert {\mathbf{\boldsymbol{x}}_1 - \mathbf{\boldsymbol{x}}_2} \rVert}. \]
Then
There exists a unique maximal solution \(\mathbf{\boldsymbol{x}}_* \colon (T_-, T_+) \to \mathbf R^n\);
Blow-up: If \(T_+\) is finite, then \(\lim_{t \to T_+} {\lVert {\mathbf{\boldsymbol{x}}(t)} \rVert} = \infty\), and similarly for \(T_-\).
\(~\)
From now on, \(t_0 = 0\) for simplicity.
Toy problem
Toy problem
For numerical illustration, we consider the equation \[ \left\{ \begin{aligned} x'(t) &= \alpha \bigl(x(t) - \sin(t) \bigr) + \cos(t) \\ x(0) &= x_0. \end{aligned} \right. \]
The exact solution is given by
\[x(t) = \sin(t) + e^{\alpha t} x_0.\]
Solutions converge if \(\alpha < 0\), and diverge if \(\alpha > 0\).
Toy problem
Toy problem
For numerical illustration, we consider the equation \[ \left\{ \begin{aligned} x'(t) &= \alpha \bigl(x(t) - \sin(t) \bigr) + \cos(t) \\ x(0) &= x_0. \end{aligned} \right. \]
The exact solution is given by
\[x(t) = \sin(t) + e^{\alpha t} x_0.\]
Solutions converge if \(\alpha < 0\), and diverge if \(\alpha > 0\).
Toy problem
Toy problem
For numerical illustration, we consider the equation \[ \left\{ \begin{aligned} x'(t) &= \alpha \bigl(x(t) - \sin(t) \bigr) + \cos(t) \\ x(0) &= x_0. \end{aligned} \right. \]
The exact solution is given by
\[x(t) = \sin(t) + e^{\alpha t} x_0.\]
Solutions converge if \(\alpha < 0\), and diverge if \(\alpha > 0\).
The forward and backward Euler methods
The forward Euler method
The method
The forward Euler method is the iteration \[ \mathbf{\boldsymbol{x}}_{n+1} = \mathbf{\boldsymbol{x}}_{n} + \Delta \mathbf{\boldsymbol{f}}(t_n, \mathbf{\boldsymbol{x}}_n), \] with \(\Delta\) the time step.
Remark: link with numerical quadrature
When \(\mathbf{\boldsymbol{f}}(t, \mathbf{\boldsymbol{x}}) = \mathbf{\boldsymbol{f}}(t)\), the method performs the approximation \[ \int_0^{n \Delta} \mathbf{\boldsymbol{f}}(t) \, \mathrm dt \approx \sum_{i=0}^{n-1} \Delta \mathbf{\boldsymbol{f}}(i \Delta). \] \({\color{green} \rightarrow}\) This is the left endpoint rule for numerical integration.
\({\color{green} \rightarrow}\) The error scales as \(\mathcal O(\Delta)\) in this case.
Question: is this true in general?
The forward Euler method: geometric interpretation
Example problem
We consider the toy problem from the beginning: \[ \left\{ \begin{aligned} x'(t) &= \alpha \bigl(x(t) - \sin(t) \bigr) + \cos(t) \\ x(0) &= 0, \end{aligned} \right. \] with \(\alpha = 0.15\).
The forward Euler method: illustration of the convergence
Estimating the rate of convergence numerically
Estimation of the order of convergence
Question: how to estimate the order \(p\)? \[ \max_n |x(t_n) - x_n| = \mathcal O(\Delta^{\color{blue} p}) \]
Answer: assume that \[ u(Δ) = C \Delta^{\color{blue} p}. \] Then \[ \log\bigl(u(\Delta)\bigr) = \log C + {\color{blue} p} \log(\Delta) \] In other words, \(\log(u)\) is an affine function of \(\log(\Delta)\)!
\({\color{green} \rightarrow}\) Enables calculating \(p\) using linear fit.
The forward Euler method: error analysis
Convergence theorem
Assume that
there exists a unique solution \(\mathbf{\boldsymbol{x}}(t) \colon [0, T] \to \mathbf R^n\), in \(C^2\) with \[ M := \sup_{t \in [0,T]} {\lVert {\mathbf{\boldsymbol{x}}''(t)} \rVert}. \]
the function \(\mathbf{\boldsymbol{f}}(t, \mathbf{\boldsymbol{x}})\) is globally Lipschitz in \(\mathbf{\boldsymbol{x}}\) with constant \(L\).
Then the following error estimate holds: \[\begin{equation} \forall n, \quad {\lVert {\mathbf{\boldsymbol{x}}(t_n) - \mathbf{\boldsymbol{x}}_n} \rVert} \leqslant \frac{{\color{blue} \Delta} M}{2} \left( \frac{\mathop{\mathrm{e}}^{L t_n} - 1}{L} \right), \end{equation}\]
Proof (dimension 1 for simplicity)
Let \(e_n = x(t_n) - x_n\). By Taylor, there is \(\tau_n \in [t_n, t_{n+1}]\) such that \[ e_n = e_{n-1} + \Delta \Bigl( f\bigl(t_{n-1}, x(t_{n-1})\bigr) - f\bigl(t_{n-1}, x_{n-1}\bigr) \Bigr) + \frac{\Delta^2}{2} x''(\tau_n). \]
By Lipschitz continuity, \[ {\lvert {e_{n}} \rvert} \leqslant(1 + \Delta L) {\lvert {e_{n-1}} \rvert} + \frac{M \Delta^2}{2}. \]
Iterating this inequality, we deduce \[ \begin{aligned} |e_{n}| &\leqslant(1 + \Delta L) \Bigl( (1 + \Delta L) {\lvert {e_{n-2}} \rvert} + \frac{M \Delta^2}{2} \Bigr) + \frac{M \Delta^2}{2}\\ &\leqslant\dotsc \leqslant(1 + \Delta L)^{n} {\lvert {e_{0}} \rvert} + \sum_{i=1}^{n} (1 + \Delta L)^{n-i} \frac{M \Delta^2}{2}. \end{aligned} \]
Using the formula for geometric series, we have \[ {\lvert {e_{n}} \rvert} \leqslant(1 + \Delta L)^{n} {\lvert {e_{0}} \rvert} + \frac{(1+\Delta L)^{n} - 1}{\Delta L} \left( \frac{\Delta^2 M}{2} \right). \]
The forward Euler method: numerical instability
Definition: numerical instability
A method is numerically unstable when the numerical solution diverges but the exact solution does not.
Example problem
Consider the toy problem from the beginning \[ \left\{ \begin{aligned} x'(t) &= f(t, x) := \alpha \bigl(x(t) - \sin(t) \bigr) + \cos(t) \\ x(0) &= x_0. \end{aligned} \right. \]
For \(\alpha < 0\), exact solutions converge with time scale \(|\alpha|^{-1}\).
For \(\alpha \ll -1\), presence of two widely separated time scales.
\({\color{green} \rightarrow}\) short scale \(|\alpha|^{-1}\), and long scale \(2\pi\) (period of \(\sin\))
\({\color{green} \rightarrow}\) the equation is called stiff.
We observe numerical instability if \(|1 + \Delta \alpha| > 1\). Why?
\({\color{green} \rightarrow}\) in general, numerical stability depends on shortest scale.
\(\alpha = -10\), \(\Delta = .18\)
\(~\)
The forward Euler method: numerical instability
Definition: numerical instability
A method is numerically unstable when the numerical solution diverges but the exact solution does not.
Example problem
Consider the toy problem from the beginning \[ \left\{ \begin{aligned} x'(t) &= f(t, x) := \alpha \bigl(x(t) - \sin(t) \bigr) + \cos(t) \\ x(0) &= x_0. \end{aligned} \right. \]
For \(\alpha < 0\), exact solutions converge with time scale \(|\alpha|^{-1}\).
For \(\alpha \ll -1\), presence of two widely separated time scales.
\({\color{green} \rightarrow}\) short scale \(|\alpha|^{-1}\), and long scale \(2\pi\) (period of \(\sin\))
\({\color{green} \rightarrow}\) the equation is called stiff.
We observe numerical instability if \(|1 + \Delta \alpha| > 1\). Why?
\({\color{green} \rightarrow}\) in general, numerical stability depends on shortest scale.
\(\alpha = -10\), \(\Delta = .22\)
\(~\)
Absolute stability
Definition: absolute stability
Consider the model equation \[ \left\{ \begin{aligned} x'(t) &= \lambda x(t), \\ x(0) &= 1. \end{aligned} \right. \] A numerical scheme is called absolutely stable if \[ \lvert x_n \rvert \xrightarrow[n \to \infty]{} 0. \] The absolute stability region is by definition \[ \mathcal A := \{ z \in \mathbf C: \text{${\lvert {x_n} \rvert} \to 0$ when $\Delta \lambda = z$} \} \subset \mathbf C. \]
Absolute stability of the forward Euler method
The forward Euler iteration for the model equation reads \[ x_{n+1} = x_n + \Delta \lambda x_n = (1 + \Delta \lambda) x_n. \] The scheme is absolutely stable iff \(|1 + \Delta \lambda| < 1\).
xlims, ylims = (-5, 1), (-2, 2)
x = LinRange(xlims..., 200)
y = LinRange(ylims..., 200)
stable = @. clamp(abs(1 + x' + im*y), 0, 2)
Plots.contourf(x, y, stable, c=cgrad(:Pastel1_3, rev=true),
aspect_ratio=:equal, levels=[0., 1., 2.],
xlabel=L"\Re(\Delta \lambda)",
ylabel=L"\Im(\Delta \lambda)",
xlims=xlims, ylims=ylims, bg=:transparent)
Plots.contour!(x, y, stable, color=:red,
levels=[0., 1., 2.], colorbar=:none)
Plots.vline!([0], color=:gray)
Plots.hline!([0], color=:gray)Absolute stability: illustration for model equation
Example problem
Consider the equation \[ \left\{ \begin{aligned} x'(t) &= - x(t), \\ x(0) &= 1. \end{aligned} \right. \] Absolute stability of forward Euler if \(\Delta < 2\).
\(\Delta = \frac{1}{2}\)
\(~\)
Absolute stability: illustration for model equation
Example problem
Consider the equation \[ \left\{ \begin{aligned} x'(t) &= - x(t), \\ x(0) &= 1. \end{aligned} \right. \] Absolute stability of forward Euler if \(\Delta < 2\).
\(\Delta = 2\)
\(~\)
Absolute stability: illustration for model equation
Example problem
Consider the equation \[ \left\{ \begin{aligned} x'(t) &= - x(t), \\ x(0) &= 1. \end{aligned} \right. \] Absolute stability of forward Euler if \(\Delta < 2\).
\(\Delta = 2.5\)
\(~\)
The main application of absolute stability: PDEs
Example: heat equation
Space discretization of the heat equation \[ \left \{ \begin{aligned} &\partial_t u = \partial_x^2 u, \\ &u(t, 0) = u(t, 1) = 0, \end{aligned} \right. \] by the finite difference method leads to the differential equation \[ \mathbf{\boldsymbol{u}}' := \mathsf A \mathbf{\boldsymbol{u}}, \qquad \mathsf A := \frac{1}{h^2} \begin{pmatrix} -2 & 1 \\ 1 & -2 & 1 \\ & 1 & -2 & 1 \\ & & \ddots & \ddots & \ddots & \\ & & & 1 & -2 & 1 \\ & & & & 1 & -2 \\ \end{pmatrix}, \] where \(h\) is the spatial step. The forward Euler method for this system reads \[ \mathbf{\boldsymbol{u}}_{n+1} = \mathbf{\boldsymbol{u}}_{n} + \Delta \mathsf A \mathbf{\boldsymbol{u}}_{n} = (\mathsf I + \Delta \mathsf A) \mathbf{\boldsymbol{u}}_{n}. \] \({\color{green} \rightarrow}\) numerical stability if \(\Delta \lambda \in \mathcal A\), i.e. \(|1 + \Delta \lambda|< 1\), for any eigenvalue \(\lambda\) of \(\mathsf A\).
\({\color{green} \rightarrow}\) since \(\max_i |\lambda_i| \propto \frac{1}{h^2}\), stability requires \(\Delta = \mathcal O(h^2)\)…
The backward Euler method
The method
The backward Euler method is the iteration \[ \mathbf{\boldsymbol{x}}_{n+1} = \mathbf{\boldsymbol{x}}_{n} + \Delta \mathbf{\boldsymbol{f}}({\color{red}t_{n+1}}, {\color{red}\mathbf{\boldsymbol{x}}_{n+1}}), \] with \(\Delta\) the time step.
\({\color{green} \rightarrow}\) Nonlinear equation to solve at each step!
\({\color{green} \rightarrow}\) A fixed point iteration is the simplest approach to solve this equation: \[ {\color{magenta}\mathbf{\boldsymbol{x}}_{n+1}^{(k+1)}} = \mathbf{\boldsymbol{x}}_{n} + \Delta \mathbf{\boldsymbol{f}}(t_n, {\color{magenta}\mathbf{\boldsymbol{x}}_{n+1}^{(k)}}). \]
Remark: link with numerical quadrature
When \(\mathbf{\boldsymbol{f}}(t, \mathbf{\boldsymbol{x}}) = \mathbf{\boldsymbol{f}}(t)\), the method performs the approximation \[ \int_0^{n \Delta} \mathbf{\boldsymbol{f}}(t) \, \mathrm dt \approx \sum_{i=1}^{n} \Delta \mathbf{\boldsymbol{f}}(i \Delta). \] \({\color{green} \rightarrow}\) This is the right endpoint rule for numerical integration.
\({\color{green} \rightarrow}\) The error scales as \(\mathcal O(\Delta)\) in this case.
Proposition absolute stability
The stability region of the backward Euler method is \[ \mathcal A = \left\{ z \in \mathbf C: \frac{1}{|1 - \Delta \lambda|} < 1 \right\} \]
\({\color{green} \rightarrow}\) Stable for all \(\Re(\Delta \lambda) < 0\), “unconditional stability”
Absolute stability: illustration for model equation
Example problem
Consider the equation \[ \left\{ \begin{aligned} x'(t) &= - x(t), \\ x(0) &= 1. \end{aligned} \right. \] Absolute stability for all \(\Delta > 0\)!
\(\Delta = 2.5\)
\(~\)
Higher-order one-step methods
Explicit Taylor methods
Motivating idea
By differentiating the equation \[ \mathbf{\boldsymbol{x}}'(t) = \mathbf{\boldsymbol{f}}\bigl(t, \mathbf{\boldsymbol{x}}(t)\bigr), \] we have access to higher-order derivatives of \(\mathbf{\boldsymbol{x}}\), e.g. \[ \begin{aligned} \mathbf{\boldsymbol{x}}'' &= \mathbf{\boldsymbol{f}}_2\bigl(t, \mathbf{\boldsymbol{x}}(t)\bigr), \qquad \mathbf{\boldsymbol{f}}_2(t, \mathbf{\boldsymbol{x}}) := (\partial_t + \mathbf{\boldsymbol{f}}(t, \mathbf{\boldsymbol{x}}) \cdot \nabla_{\mathbf{\boldsymbol{x}}}) \mathbf{\boldsymbol{f}}(t, \mathbf{\boldsymbol{x}}) \\ \mathbf{\boldsymbol{x}}''' &= \mathbf{\boldsymbol{f}}_3\bigl(t, \mathbf{\boldsymbol{x}}(t)\bigr), \qquad \mathbf{\boldsymbol{f}}_3(t, \mathbf{\boldsymbol{x}}) := (\partial_t + \mathbf{\boldsymbol{f}}(t, \mathbf{\boldsymbol{x}}) \cdot \nabla_{\mathbf{\boldsymbol{x}}}) \mathbf{\boldsymbol{f}}_2(t, \mathbf{\boldsymbol{x}}), \\ &\dots \end{aligned} \] This motivates the so-called Taylor methods: \[ \mathbf{\boldsymbol{x}}_{n+1} = \mathbf{\boldsymbol{x}}_n + \Delta \mathbf{\boldsymbol{f}}(t_n, \mathbf{\boldsymbol{x}}_n) + \frac{\Delta^2}{2} \mathbf{\boldsymbol{f}}_2(t_n, \mathbf{\boldsymbol{x}}_n) + \dotsc + \frac{\Delta^p}{p!} \mathbf{\boldsymbol{f}}_p(t_n, \mathbf{\boldsymbol{x}}_n). \]
Remarks
Note that for \(p = 1\), this is the forward Euler method;
The error scales as \(\mathcal O(\Delta^p)\) under appropriate assumptions.
By construction, exact when \(\mathbf{\boldsymbol{x}}(t)\) is polynomial of degree \(p\).
Proposition: absolute stability
The absolute stability region of the Taylor method is \[ \mathcal A = \left\{ z \in \mathbf C: \left|1 + z + \dotsc + \frac{z^p}{p!} \right| < 1 \right\} \]
Runge-Kutta methods
Definition: explicit Runge-Kutta method
An explicit Runge–Kutta method with \(p\) stages is of the form \[ \begin{aligned} \mathbf{\boldsymbol{x}}_{n+1} &= \mathbf{\boldsymbol{x}}_n + \Delta \sum_{i=1}^p b_i \mathbf{\boldsymbol{k}}_i \\ \mathbf{\boldsymbol{k}}_1 &= \mathbf{\boldsymbol{f}}(t_n, \mathbf{\boldsymbol{x}}_n), \\ \mathbf{\boldsymbol{k}}_2 &= \mathbf{\boldsymbol{f}}\bigl(t_n + c_2 \Delta, \mathbf{\boldsymbol{x}}_n+\Delta(a_{21} \mathbf{\boldsymbol{k}}_1)\bigr), \\ \mathbf{\boldsymbol{k}}_3 &= \mathbf{\boldsymbol{f}}\bigl(t_n + c_3 \Delta, \mathbf{\boldsymbol{x}}_n+\Delta(a_{31} \mathbf{\boldsymbol{k}}_1 + a_{32} \mathbf{\boldsymbol{k}}_2)\bigr), \\ &\;\;\vdots \\ \mathbf{\boldsymbol{k}}_p &= \mathbf{\boldsymbol{f}}\left(t_n + c_p \Delta, \mathbf{\boldsymbol{x}}_n + \Delta \sum_{j = 1}^{p-1} a_{pj} \mathbf{\boldsymbol{k}}_j\right), \end{aligned} \]
\({\color{green} \rightarrow}\) Derivatives of \(\mathbf{\boldsymbol{f}}\) not necessary!
Remarks
Construction tedious and not discussed here;
Most widely used in applications;
Error scales as \(\mathcal O(\Delta^p)\) for appropriate coefficients;
Example: Heun’s method
Heun’s method is a Runge-Kutta method of order 2 and reads \[ \mathbf{\boldsymbol{x}}_{n+1} = \mathbf{\boldsymbol{x}}_n + \frac{\Delta}{2}\mathbf{\boldsymbol{f}}(t_n, \mathbf{\boldsymbol{x}}_n) + \frac{\Delta}{2} \mathbf{\boldsymbol{f}} \bigl( t_n + \Delta, \mathbf{\boldsymbol{x}}_n + \Delta \mathbf{\boldsymbol{f}}(t_n, \mathbf{\boldsymbol{x}}_n) \bigr). \]
Proposition: same stability region as Taylor’s method
The absolute stability region of a Runge-Kutta method of order \(p\) is \[ \mathcal A = \left\{ z \in \mathbf C: \left|1 + z + \dotsc + \frac{z^p}{p!} \right| < 1 \right\} \]
Further readings
Some of the topics we did not cover include
Multistep methods;
Implicit Runge-Kutta methods;
Adaptive time-stepping;
Symplectic schemes for Hamiltonian systems;
…
