Notes - Numerical Analysis HT24, Initial value problems
Overview of this part of the course
In this part of the course, we learn about ODE initial value problems:
\[\pmb y' (x) = \pmb f(x, \pmb y), \quad \pmb y(x_0) = \pmb y\]Most ODEs don’t have analytic solutions, so we need numerical methods to solve them. We consider:
- One-step methods (like Euler’s method, we update our approximation in a direction informed by the one evaluation of the derivative, i.e. function $\pmb f$)
- Runge-Kutta methods (e.g. RK4, we also move in a direction informed by the value of the derivative, but evaluate $\pmb f$ at lots of places between each $x _ i$ and $x _ {i+1}$)
- Multi-step methods (e.g. Adams–Bashforth, we also move in a direction informed by the value of the derivative, but evaluate $\pmb f$ at fixed grid points between each $x _ i$ and $x _ {i+1}$)
Each method has a step size $h$, which control how small our steps are between each approximation. We evaluate each numerical method by its:
- Convergence: Does the method converge to the exact solution of the ODE as $h$ tends towards $0$? If so, how fast? What is the “order of accuracy” (the exponent of $h$ in the expression $e = O(h^p)$ where $e$ is the error at the end of the interval)?
- Stability: How small does $h$ need to be before it starts looking like the exact solution? When does the solution explode off to infinity when it really should converge to $0$?
In summary, the analysis of each method under these benchmarks roughly shows:
- For one-step methods:
- Convergence: any consistent one-step method is convergent, where consistency roughly means that the local error at each approximation can be made arbitrarily small by shrinking the step size. If the consistency order is $p$, then the order of convergence is also $p$.
- Stability: since one-step methods can be seen as a type of Runge-Kutta method, this is covered by the analysis of Runge-Kutta methods
- For Runge-Kutta methods:
- Convergence: any consistent Runge-Kutta method is convergent where consistency is defined in the same way as for one-step methods, and the results are very similar to one-step methods.
- Stability: one benchmark of performance is whether the method does a “good job” on the ODE $y’ = zy$, $y(0) = 1$, which has solution $y(x) = e^{zx}$, and where “good job” means that if $\text{Re }z < 0$, then $y _ n \to 0$ as $n \to \infty$ (where $y _ n$ is the approximation at $y(x _ n)$). To evaluate this, we relate $y _ n$ to $S(z)$, the “stability function” of the method, defined as $S(z) = \pmb \Psi(0, 1, 1, y \mapsto yz)$. By unrolling some recurrences, it turns out that a sufficient condition to do a “good job” is that $\mathbb C^- \subseteq S _ \Psi$ (so called $A$-stability), where $\mathbb C^-$ is the open left half-plane and $S _ \Psi$ is the stability region defined as the set of $z$ for which $ \vert S(z) \vert < 1$. It also is useful to have $L$-stability, which is the condition that $\lim _ {\text{Re}(z) \to -\infty} \vert S(z) \vert = 0$, and roughly corresponds to saying that the method converges faster and faster to $0$ the more negative $z$ is in $e^{zx}$.
- For multi-step methods:
- Convergence: Convergence turns out to be more complicated then for one-step and Runge-Kutta methods. We have a slightly different notion of what it means to be consistent, and there is a big theorem that states that consistency and the so-called “root condition” (the associated characteristic polynomial of the method has all roots in the closed unit disc, and all roots on the boundary are simply) are necessary and sufficient for convergence.
- Stability: Analogously to Runge-Kutta methods, we also consider the performance of the methods on the Dahlquist test equation for different values of $z$. By again relating what the method does with step size $1$ and squinting hard enough to see a recurrence relation, it turns out that the analog of the stability function is the stability polynomial $\pi(x; z)$. This motivates the definition of the stability region of a multistep method as the set of all points $z \in \mathbb C$ such that if $x$ is a root of $\pi$ (i.e. $\pi(x, z) = 0$), then it must lie in the closed unit disc, and if it is not a simple root, then $ \vert x \vert < 1$. Then there are identical definitions of $A$-stability and $L$-stability.
In more detail:
- One-step methods:
- Some function $\pmb \Psi(s, \pmb y, h, \pmb f)$ which computes an approximation $\pmb y(s + h)$ to the initial value problem specified by $\pmb y’ = f$ and $\pmb y(s) = \pmb y$ (the notation is a bit sloppy, since $\pmb y$ is both the name of the function and the name of the intial point)
- Included among them is Euler’s method, the midpoint method, etc
- To analyse convergence, we use the idea of “consistency”
- A one-step method is consistent if $\pmb \Psi(s, \pmb y, 0, \pmb f) = \pmb y$ and $\frac{\text d}{\text d h} \pmb \Psi(s, \pmb y, h, \pmb f) \Big \vert _ {h = 0} = \pmb f(s, \pmb y)$.
- The consistency error is defined by
- A method is consistent iff the consistency error tends towards $0$ as $h$ tends towards $0$
- Consistency order is defined analogously to the order of accuracy, it is the exponent of $h$ in the expression $\pmb \tau(s, \pmb y, h, \pmb f) = O(h^p)$
- *Convergence*: There are theorems relating the global error of the approximation and the consistency error, and it turns out that as long as the consistency error tends to $0$ as $h$ tends to $0$ (i.e. the method is consistent), then the method converges
- *Stability*: Since one-step methods are a special case of Runge-Kutta methods, this is analysed in a later section. - Runge-Kutta methods
- Convergent one-step methods make approximations via $\pmb y _ {n+1} = \pmb y _ n + h\pmb f(\cdots)$ (this is not the definition of one-step methods, but follows from a theorem relating consistency to continuous increment functions)
- Runge-Kutta methods are more sophisticated, they instead make approximations via $\pmb \Psi(s, \pmb y, h, \pmb f) = \pmb y + h \sum^s _ {i = 1} b _ i \pmb k _ i$ where $\pmb k _ i$ are the "stages" $\pmb k _ i := \pmb f(x + c _ i h, \pmb y + h \sum^s _ {j = 1} a _ {ij} \pmb k _ j)$ and $c _ i = \sum^s _ {j = 1} a _ {ij}$, i.e. we use multiple evaluations of the function at different points between $x _ n$ and $x _ {n+1}$ in order to get a better sense of which direction the differential equation is headed in
- Consistency error still makes sense for Runge-Kutta methods
- Various results about the consistency of Rugne-Kutta methods:
- Consistent (i.e. at least consistency order 1) iff $\sum^s _ {i = 1} b _ i = 1$
- At least consistency order 2 if above and $\sum^s _ {i = 1} b _ i c _ i = \frac 1 2$
- At least consistency order 3 if all of the above and $\sum _ {i=1}^{s} b _ i \sum _ {j=1}^{s} a _ {ij} c _ j = \frac{1}{6}$ and $\sum _ {i=1}^{s} b _ i c _ i^2 = \frac{1}{3}$
- Lots and lots of conditions for higher order consistency, e.g. in order to be order 20, there are $20247374$ conditions
- The consistency order $p$ is bounded in general by twice the number of stages $s$
- If the method is explicit, then the consistency order $p$ is bounded by the number of stages $s$.
- To analyse stability, we consider the Dahlquist test equation $y' = zy$, $y(0) = 1$ and $\text{Re}(z) < 0$ which has solution $y(x) = \exp(xz)$ and $y(x) \to 0$ as $x \to \infty$. We want Runge-Kutta methods to respect this, and have the approximation $y _ n \to 0$ as $n \to \infty$.
- We then associate with each Runge-Kutta method $\pmb \Psi$ a stability function $S : \mathbb C \to \mathbb C$, which, computes the approximation that $\pmb \Psi$ would make with step size $1$ on the the Dahlquist test equation, i.e. $S(z) = \pmb \Psi(0, 1, 1, y \mapsto yz)$
- Then we show that if we use that Runge-Kutta method to get an approximate solution $y _ n$ to the Dahlquist test equation with value $z$, $y _ n = S(zh)^n$. So it then follows that for the approximate solution given by the Runge-Kutta to converge to $0$, we need $S(zh) < 1$.
- This motivates the definition of the stability region $S _ {\pmb \Psi}$ as $\\{z \in \mathbb C \mid S(z) < 1\\}$.
- We care about when $\text{Re}(z) < 0$, since this is when the corresponding Dahlquist test equation converges to $0$
- A method is called $A$-stable if $\mathbb C^{-} = \\{z \in \mathbb C \mid \text{Re}(z) < 0\\}$
- $A$-stability ensures that the solution $y _ n$ eventually tends to $0$, but it could be that this happens very slowly. If $z = -1000$, the exact solution of the Dahlquist test equation would converge to $0$ rapidly, but an $A$-stable Runge-Kutta method might converge very slowly
- This motivates the definition of $L$-stability. If an $A$-stable method also satisfies $\lim _ {\text{Re}(z) \to -\infty} \vert S(z) \vert = 0$, then it is said to be $L$-stable - Multi-step methods
- Runge-Kutta methods assume that we can evaluate $\pmb f$ anywhere, this might not be true in practice
- Multi-step methods only use evaluations of $\pmb f$ at grid points
- They compute $\pmb y _ {n+k}$ via the approximation $\sum^k _ {j = 0} \alpha _ j \pmb y _ {n + j} = h \sum^k _ {j = 0} \beta _ j \pmb f(x _ {n + j}, \pmb y _ {n + j})$
- Summarise each method using a polynomial: $\rho(z) = \sum^k _ {j = 0} \alpha _ j z^j$ and $\sigma(z) = \sum^k _ {j = 0} \beta _ j z^j$.
- Quadrature rules naturally give rise to multistep methods: this can be seen by writing out the approximation used in rectangle rules and then applying the quadrature rule to $\pmb f$ in the integral
- *Convergence*: Can't define consistency order of a multi-step method as there's no simple function that gives "the next approximation". Instead we define the consistency error as $\pmb \tau(h) = \frac{\sum^k _ {j = 0} \alpha _ j \pmb y(x _ j) - h \sum^k _ {j = 0} \beta _ j \pmb y'(x _ j)}{h \sum^k _ {j = 0}\beta _ j}$ (I am not sure how to motivate this definition, though it reduces to the definition for one-step methods when $k = 1$ and $\alpha _ 1 = 1$)
- Taylor expansion of $\pmb y$ gives us the result that a multi-step method is consistent iff $\rho(1) = 0$ and $\rho'(1) = \sigma(1) \ne 0$
- However, unlike in the case with one-step and Runge-Kutta methods, there isn't as straightforward a relationship between consistency and convergence; it's possible to come up with consistent methods that don't converge, one way to see this is to consider when $h$ is very small, the definition effectively behaves like a recurrence relation, and the exponents involved in the exact solution depend on the roots of the first characteristic polynomial $\rho(z)$.
- The root condition is: all roots of $\rho(z)$ lie on the closed unit disc and all roots that lie on the boundary are simple. This means that the recurrence won't blow up (why the condition about the roots on the boundary being simple? consider the Jordan blocks if the root is not simple, these will blow up)
- There is then a big theorem called Dahlquist's equivalence theorem which states that consistency and the root condition (also called zero-stability) is enough for convergence, and that the starting values are consistent, then the global error order (and in fact the error over the whole interval) is the same as the consistency order.
- *Stability*: Use the same approach as Runge-Kutta methods, where we analyse what the solution looks like on the Dahlquist test equation. It turns out that the $n$-th iterate looks like $\sum^k _ {j = 0} \alpha _ j y _ {n+j} = \sum^k _ {j = 0} \beta _ j z y _ {n+j}$, or equivalently $\sum^k _ {j = 0} (\alpha _ j - \beta _ j z) y _ {n+j} = 0$ . Then solving the recurrence gives $y _ n = p _ 1(n) r^n _ 1 + \cdots + p _ \ell(n) r^n _ \ell$ where the $r _ j$ are the roots of the polynomial $\pi(x; z) = \sum^k _ {j = 0} (\alpha _ j - \beta _ j z) x^j = 0$ and $p _ j$ are polynomials of degree $m _ j$ and $m _ j$ is the multiplicity of $r _ j$.
- Then:
- if $\pi(x; z)$ has a zero outside the unit disc, $y _ n$ grows exponentially
- if any $r _ j$ is on the unit disc and is not simple, then $y _ n \sim n^{m _ j -1}$
- otherwise (all roots of $\pi(x; z)$ are in the unit disc), $y _ n$ converges to $0$
- $\pi(x; z)$ is called the stability polynomial, can also be written $\pi(x; z) = \rho(x) - z\sigma(x)$.
- This motivates the definition of the stability region of a multistep method as
Flashcards
Define an first-order ordinary differential equation (ODE) and a subsequent initial value problem (IVP).
Let
- $I = [x _ 0, X] \subset \mathbb R$ be an interval
- $D \subset \mathbb R^d$ open subset
An ODE is an equation of the form
\[\pmb y' (x) = \pmb f(x, \pmb y)\]where $\pmb f : I \times D \to \mathbb R^d$. An initial value problem adds the extra condition that
\[\pmb y(x_0) = \pmb y_0\]The goal is to compute $\pmb y(x)$ for different values of $x$.
Suppose:
- $I = [x _ 0, X] \subset \mathbb R$ be an interval
- $D \subset \mathbb R^d$ open subset
- $\pmb y’ (x) = \pmb f(x, \pmb y)$
- $\pmb y (x _ 0) = \pmb y _ 0$
Then can you give a closed form for $\pmb y(x)$ (though not effective for calculation)?
\[\pmb y(x) = \pmb y(x_0) + \int^x_{x_0} \pmb f(t, \pmb y(t)) \text d t\]
Suppose:
- $I = [x _ 0, X] \subset \mathbb R$ be an interval
- $D \subset \mathbb R^d$ open subset
- $\pmb y’ (x) = \pmb f(x, \pmb y)$
- $\pmb y (x _ 0) = \pmb y _ 0$
- $\pmb f$ is $L$-Lipschitz, i.e. $ \vert \vert \pmb f(x, \pmb y) - \pmb f(x, \pmb z) \vert \vert \le L \vert \vert \pmb y - \pmb z \vert \vert $ $\forall (x, \pmb y)$ and $(x, \pmb z)$.
Then can you define what it means for solutions to be “stable”?
If
- $\pmb y : [x _ 0, X] \to D$ solves the IVP with $\pmb f(x _ 0) = \pmb y _ 0$
- $\hat{\pmb y} : [x _ 0, X] \to D$ solves the IVP with $\pmb f(x _ 0) = \hat{\pmb y} _ 0$
Then $\forall x \in [x _ 0, X]$:
\[||\pmb y (x) - \hat{\pmb y}(x)|| \le e^{L(X - x_0)} ||\pmb y_0 - \hat{\pmb y_0}||\]Suppose:
- $I = [x _ 0, X] \subset \mathbb R$ be an interval
- $D \subset \mathbb R^d$ open subset
- $\pmb y’ (x) = \pmb f(x, \pmb y)$
What does it mean for the ODE to be autonomous, and how can you reformulate any ODE as an autonomous one?
Autonomous: $\pmb f$ does not depend on $x$.
We can just construct $\hat{\pmb y} = [x, \pmb y]^\top$ and then consider
\[\hat{\pmb y}' = \pmb f(\hat{\pmb y})\]