# Stability of timestepping schemes

You already saw in CMIII that for some problems if we make the timestep too large in the forward Euler method, we don't get a valid solution. Let's remind ourselves with an example.

## Stiff problems

Let's consider a seemingly benign linear system of equations (this example is taken from Iserles' book, section 4.1)

$$
\begin{bmatrix}
\dot u_0\\
\dot u_1\\
\end{bmatrix}
= \overbrace{\begin{bmatrix}
-100 & 1\\
0 & -1/10\\
\end{bmatrix}}^{\Lambda}
\begin{bmatrix}
u_0\\
u_1\\
\end{bmatrix}.
$$

Factorising the matrix $\Lambda = V D V^{-1}$, it is easy to show that the exact solution at time $t$ obeys

$$
u(t) = e^{-100 t} x_1 + e^{-t/10} x_2
$$

where $x_1$ and $x_2$ depend on the initial conditions, but not $t$.

For large $t$, the behaviour of the system is controlled by the slowest decaying of the two exponentials. That is, for all practical purposes, we can write:

$$
u(t) \approx e^{-t/10} x_2, t > 0.
$$

Let's try out explicit Euler for this problem.


```python
%matplotlib widget
import numpy
from matplotlib import pyplot
pyplot.style.use("ggplot")

from scipy.linalg import expm

class linear(object):
    def __init__(self, A):
        self.A = numpy.asarray(A)
    def f(self, t, u):
        return self.A @ u
    __call__ = f
    def explicit(self, t, u):
        return numpy.zeros_like(u)
    def J(self, t, u):
        return self.A
    def exact(self, t, u_0):
        t = numpy.array(t, ndmin=1)
        return [numpy.real_if_close(expm(self.A*s) @ u_0) for s in t]
    
test = linear([[-100, 1], [0, -1/10]])

u_0 = numpy.asarray([1, 1])

def ode_euler(f, u_0, h=0.1, T=1):
    u = numpy.array(u_0)
    t = 0
    thist = [t]
    uhist = [u_0]
    while t < T:
        h = min(h, T - t)
        u = u + h * f(t, u)
        t += h
        thist.append(t)
        uhist.append(u.copy())
    return numpy.asarray(thist), numpy.asarray(uhist)

thist, uhist = ode_euler(test, u_0, h=1/100, T=1)
pyplot.figure()

pyplot.semilogy(thist, uhist, "o", linestyle="solid", label='Forward Euler')
pyplot.semilogy(thist, test.exact(thist, u_0), label='exact');
pyplot.legend();
```

### Observations

When the timestep is very small, we get a reasonably accurate answer. However, when the timestep gets big, suddently the numerical solution blows up.

This is due to the small _stability domain_ of the explicit Euler method. As we will see below, the problem is only stable in a small region in the complex plane.

This constraint forces us to _artificially_ constrain the step size to avoid instabilities in the numerical solution.

We will informally call these problems, where the timestep is small for stability, rather than accuracy, reasons *stiff*

### Questions

1. Show that the explicit Euler solution after $n$ steps of size $h$ can be written as
$$
u_n = (1 - 100h)^n x_1 + (1 - \frac{1}{10}h)^n x_2
$$
2. Use this to determine the appropriate maximum step size for the integration.
3. Confirm your analysis with some experiments

## An implicit method

Now let's consider the same problem, only this time we'll use backward Euler, an implicit method.


```python
def ode_beuler(eqn, u_0, h=0.1, T=1):
    A = eqn.A
    u = numpy.array(u_0)
    t = 0
    thist = [t]
    uhist = [u_0]
    while t < T:
        h = min(h, T - t)
        # u <- (I - h A)^{-1} u
        u = numpy.linalg.solve(numpy.eye(len(A)) - h*A, u)
        t += h
        thist.append(t)
        uhist.append(u.copy())
    return numpy.asarray(thist), numpy.asarray(uhist)
```


```python
thist, uhist = ode_beuler(test, u_0, h=1, T=10)
pyplot.figure()

pyplot.semilogy(thist, uhist, "o", linestyle="solid", label='Backward Euler')
pyplot.semilogy(thist, test.exact(thist, u_0), label='exact');
pyplot.legend();
```


```python

```

Now we can take large timesteps, and only need to make $h$ small if we want higher accuracy.

## Linear stability

Why is this?

To answer this, we consider the behaviour of the integrator on a linear test problem, the *Dahlquist test equation*:

$$
\dot u = \lambda u \quad \lambda \in \mathbb{C}.
$$

When taking a step of length $h$, this equation has the exact solution

$$
u(h) = u_0 e^{\overbrace{\lambda h}^z} = u_0 e^{\mathfrak{R} z}(\cos \mathfrak{I} z + i\sin\mathfrak{I} z),
$$
where we wrote $\mathbb{C} \ni z = \lambda h$.

This problem is _physically_ stable whenever the real part of $z$ (and hence $\lambda$) is less than zero: $\mathfrak{R} z \le 0$.

### Aside

There is also theory surrounding _nonlinear_ stability, but we will not cover it in this course.

## Stability regions

Let us now consider applying the timestepping schemes we've already encountered to the test equation to see how they behave.

### Explicit Euler

$$
u(h) = \underbrace{(1 + \lambda h)}_{R(z)}u_0.
$$

Repeated application of the scheme results in:

$$
u(mh) = R(z)^m u_0.
$$

This scheme is convergent (producing a bounded $u(mh)$ in the limit $m \to \infty$) only if

$$
|R(z)| \le 1.
$$

### Implicit Euler

Performing a similar calculation we obtain

$$
R(z) = \frac{1}{1 - z}.
$$

### Trapezoidal rule/implicit midpoint

$$
R(z) = \frac{1 + z/2}{1 - z/2}.
$$

### Definitions
The function $R(z)$ is called the *stability function* of the method, the set

$$
S = \{ z \in \mathbb{C} \colon |R(z)| \le 1 \}
$$

is called the *stability domain*.

## Stability plots

Visualising the stability domains, by plotting $|R(z)|$ is an insightful way of comparing methods:


```python
def plot_stability(x, y, R, label):
    pyplot.figure()
    C = pyplot.contourf(x, y, numpy.abs(R), numpy.linspace(0, 1, 10), cmap=pyplot.cm.coolwarm)
    
    pyplot.colorbar(C, ticks=numpy.linspace(0, 1, 10))
    pyplot.contour(x, y, numpy.abs(Rz), numpy.linspace(0, 1,4), colors='k')
    pyplot.xlabel("Re (z)")
    pyplot.ylabel("Im (z)")
    pyplot.title(label)
```


```python
x = numpy.linspace(-3, 3)
x, y = numpy.meshgrid(x, x)
z = x + 1j*y

Rs = [("Forward Euler", 1 + z),
      ("Backward Euler", 1/(1 - z)),
      ("Implicit Midpoint", (1 + z/2)/(1 - z/2))]

Rs = [("Heun", 1 + z * (1 + z/2))]
for label, Rz in Rs:
    plot_stability(x, y, Rz, label)
```


```python
pyplot.figure()
z = numpy.linspace(-2, 1)
pyplot.plot(z, 1 + z*(1 + z/2))
```

### Observations

While the physical problem is stable whenever $\mathfrak{R} \lambda \le 0$, the same is not true of all our methods. The explicit Euler method is only stable in a small region of the left half plane, implicit Euler is stable in the entire left half plane, and more, implicit midpoint is stable in exactly the left half plane.

A question naturally arises? Why would one ever use an explicit method? Can you think of any reasons?

## Definitions

### A-stability

A method is _A-stable_ if the stability domain

$$
S = \{z \in \mathbb{C} \colon |R(z)| \le 1\}
$$

contains the _entire_ left half plane

$$
\mathfrak{R} z \le 0.
$$

This means that we can take arbitrarily large timesteps ($h \to \infty$) without the method becoming unstable (diverging) for any problem that is physically stable. Note that this says nothing about the _accuracy_ of the method.

No explicit method is A-stable.

### Questions

1. Show that the midpoint method and implicit Euler really do contain the entire left half plane. Hint: multiply through by an appropriate choice of $1$, write $z = a + bi$ and show that $|R(z)| \le 1$ whenever $a \le 0$.


```python

```

# Higher order methods

So far, we've seen the implicit and explicit Euler methods which are only first order accurate, and the implicit midpoint method which is second order accurate. We might, especially if we have an accurate spatial discretisation, want more accuracy in our time integrator. We will now look at some methods that offer this. Certainly if the timestep size is no longer limited by stability requirements (for example when using an implicit integrator), it makes sense to take big steps, but then we need a high order method for accuracy.

In this section we will look at single-step Runge-Kutta methods (covered briefly at the very end of CMIII).

Many of these methods offer a clever way of estimating the local error accumulation at each step, and we will look at how this can be used to implement an adaptive timestepping scheme.

High-order timestepping schemes also have a higher cost of evaluation (as we shall see), and we'll discuss how to compare different schemes in a fair way using _work-precision_ diagrams.

## [Runge-Kutta methods](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods)

These are a class of _single-step_ methods. That is, we advance a solution $u_0$ from $t_0$ to $t_1$ using $u_0$ as an initial value, producing $u_1$, then when advancing from $t_1$ to $t_2$ we forget about $u_0$ and use $u_1$ as an initial value.

Recall that the single-step methods we have looked at so far (Euler and midpoint schemes) are all obtained by formally solving the ODE with integration and then approximating the integral.

$$
\dot u = f(t, u), u(t_0) = u_0
$$

and hence

$$
u(t) = u(t_0) + \int_{t_0}^t f(\tau, u) \text{d}\tau.
$$

A single step is then obtained by setting $t = t_0 + h$ and writing

$$
\int_{t_0}^{t_0 + h} f(\tau, u) \text{d}\tau = h\int_0^1 f(t_0 + h\tau, y(t_0 + h\tau))\text{d}\tau.
$$

We now _approximate_ the integral with a sum (using a quadrature scheme, but let's sweep that under the carpet for now), by picking a set of points $\{\xi_i\}_{i=1}^{N} \in [0, 1]$ and weights $\{w_i\}_{i=1}^{N} \in \mathbb{R}$

$$
\int_0^1 f(t_0 + h\tau, y(t_0 + h\tau))\text{d}\tau \approx \sum_{i=1}^{N} w_i f(t_0 + \xi_i h, y(t_0 + \xi_i h)).
$$

Giving us our update formula:
$$
u_{n+1} = u_n + h \sum_{i=1}^{N} w_i f(t_n + \xi_i h, u(t_n + \xi_i h)).
$$

For example, explicit Euler was obtained with $N=1$, $\{\xi_i\} = \{0\}$, and $\{w_i\} = \{1\}$.

Our challenge, for more general schemes, is to figure out how to approximate the (unknown) $u$s, we only know the initial value $u(t_0)$.

### Explicit Runge-Kutta methods

Write $Y_i$ for our numerical approximation to $u(t_0 + \xi_i h)$. _Explicit_ Runga-Kutta methods compute $Y_i, i = 2,\dots,N$ as a linear combination of $Y_j, j < i$.

That is, we write

$$
\begin{align}
Y_1 &= u_n\\
Y_2 &= u_n + h a_{2, 1} f(t_n, Y_1)\\
Y_3 &= u_n + h(a_{3,1} f(t_n, Y_1) + a_{3,2} f(t_n + c_2 h, Y_2)\\
&\vdots\\
Y_N &= u_n + h \sum_{i=1}^{N-1} a_{N, i} f(t_n + c_i h, Y_i)\\
u_{n+1} &= u_n + h \sum_{i=1}^{N} b_{i} f(t_n + c_i h, Y_i).
\end{align}
$$

Each of these intermediate steps is termed a *stage*, and so we just wrote down an $N$-stage Runge-Kutta method.

For an explicit method, the _RK matrix_ (or *coefficients*)

$$
A = \begin{bmatrix}
a_{1, 1} & a_{1, 2} & \dots & a_{1, N}\\
a_{2, 1} & a_{2, 2} & \dots & a_{2, N}\\
\vdots & \vdots & \ddots & \vdots\\
a_{N, 1} & a_{N, 2} & \dots & a_{N, N}\\
\end{bmatrix}
$$

is strictly lower-triangular ($a_{i, j} = 0, j \ge i$).

If we further write the _RK weights_ (or *completion weights*)

$$
b = \begin{bmatrix}
b_1 \\
b_2 \\
\vdots \\
b_N
\end{bmatrix}
$$

and the _RK nodes_ (or *abscissa*)

$$
c = \begin{bmatrix}
c_1 \\
c_2 \\
\vdots \\
c_N
\end{bmatrix}
$$

the entire scheme can be compactly represented in a _Butcher tableau_ 

$$
\begin{array}{c|c}
c & A\\
\hline
 & b^T
\end{array}
$$

(named after [John Butcher](https://en.wikipedia.org/wiki/John_C._Butcher) who developed much of the theory underpinning the analysis and design of Runge-Kutta methods).

### Determining $A$, $b$, and $c$

For low order methods, our first thought is to Taylor expand everything around $(t_n, u_n)$ and then equate terms with the Taylor expansion of the exact solution.

This is doable for $N=2$, but already becomes tremendously labourious at $N=3$. Fortunately, other people (notably Butcher) have developed graph theoretical tools to construct, and analyse the properties of, Runge-Kutta methods. Here, we will only the resulting methods.

One particular point to note is that the conditions on the entries in $A$, $b$, and $c$ do not uniquely specify the coefficients, so it is possible to come up with multiple different methods with the same number of stages.

After all that, let's write some code to do explict RK integration. The Butcher tableau is not just a useful mathematical device, it allows us to write generic code for RK integration (rather than one method for explicit Euler, one for implicit Euler, one for midpoint, etc...).


```python
class ButcherTable(object):
    def __init__(self, A, b, order=None, btilde=None):
        self.A = numpy.asarray(A)
        self.b = numpy.asarray(b)
        self.order = order
        self.btilde = numpy.asarray(btilde) if btilde else None
        self.c = self.A.sum(axis=1)

    def _repr_latex_(self):
        rows, cols = self.A.shape
        strs = ["$$", r"\left[", r"\begin{array}{c|%s}" % ("c"*cols)]
        for r in range(rows):
            row = " & ".join(map(str, [self.c[r]] + list(self.A[r, :])))
            strs.append(row)
            strs.append(r"\\")
        strs.append(r"\hline")
        
        strs.append("&")
        strs.append(" & ".join(map(str, self.b)))
        strs.append(r"\\")
        if self.btilde is not None:
            strs.append("&" + " & ".join(map(str, self.btilde)))
        strs.extend([r"\end{array}", r"\right]", "$$"])
        return "\n".join(strs)
```


```python
euler = ButcherTable([[0]], [1], order=1); euler
```


```python
beuler = ButcherTable([[1]], [1], order=1); beuler
```


```python
midpoint = ButcherTable([[1/2]], [1], order=2); midpoint
```


```python
def ode_rk_explicit(f, u_0, table, h=0.1, T=1):
    A = table.A
    b = table.b
    c = A.sum(axis=1)
    N = len(c)
    for i in range(N):
        if any(A[i, i:] > 0):
            raise ValueError("Only for explicit methods")
    t = 0
    u_0 = numpy.asarray(u_0)
    u = u_0.copy()
    hist = [(t, u.copy())]
    fY = numpy.zeros(u_0.shape + (N, ))
    while t < T:
        h = min(h, T - t)
        for i in range(N):
            Yi = u.copy()
            for j in range(i):
                Yi += h * A[i, j] * fY[:, j]
            fY[:, i] = f(t + h*c[i], Yi)
        u += h * fY @ b
        t += h
        hist.append((t, u.copy()))
    ts, us = zip(*hist)
    return numpy.asarray(ts), numpy.asarray(us)
```


```python
test = linear([[0, 1.0], [-1.0, 0.0]])

u_0 = numpy.asarray([1.0, 0])

ts, us = ode_rk_explicit(test, u_0, euler, h=0.1, T=15)
pyplot.figure()
pyplot.plot(ts, us, label="Forward Euler");
pyplot.legend();
```

We can also combine multiple explicit Euler steps for a higher order method, [Heun's method](https://en.wikipedia.org/wiki/Heun%27s_method) (also called "modified Euler" or "explicit trapezoid"). This has Butcher tableau

$$
\left[
\begin{array}{c|cc}
0 & 0 & 0\\
1 & 1 & 0\\
\hline
  & 1/2 & 1/2
\end{array}
\right]
$$


```python
heun = ButcherTable(numpy.asarray([[0, 0], [1, 0]]), numpy.asarray([1/2, 1/2]), order=2); heun
```


```python
ts, us = ode_rk_explicit(test, u_0, heun, h=0.1, T=15)
pyplot.figure()
pyplot.plot(ts, us, label="Heun");
pyplot.legend();
```

Perhaps _the_ most famous RK method is the explicit 4-th order RK4


```python
rk4 = ButcherTable(numpy.asarray([[0, 0, 0, 0],
                                  [1/2, 0, 0, 0],
                                  [0, 1/2, 0, 0],
                                  [0, 0, 1, 0]]),
                   numpy.asarray([1/6, 1/3, 1/3, 1/6]), order=4); rk4
```


```python
ts, us = ode_rk_explicit(test, u_0, rk4, h=0.1, T=15)
pyplot.figure()
pyplot.plot(ts, us, label="RK4");
pyplot.legend();
```

## Questions

1. Derive an expression for the stability function $R(z)$ for RK methods by applying one step to the Dahlquist test equation
2. Extend the `ButcherTable` class with a method to plot the method's stability region.
3. What do you see when comparing the stability of Heun and RK4 to that of explicit Euler?


```python

```

## Comparing methods

All other things being equal, a higher-order method is likely a superior choice, since it will converge faster to the exact solution with bigger timesteps. For example, let's compute the convergence rate for the explicit Euler, BS3, RK4, and DP5 methods on our linear osciliatory test problem.


```python
from functools import partial

def mms(eqn, integrator, u_0, h=0.1, T=15):
    thist, uhist = integrator(eqn, u_0, h=h, T=T)
    return numpy.linalg.norm(uhist[-1] - eqn.exact(thist[-1], u_0), numpy.inf)

eqn = linear(numpy.array([[0, 1],
                          [-1, 0]]))
u_0 = numpy.array([.75, 0])

hs = numpy.array([1/2**i for i in range(2, 10)])
pyplot.figure()
method_errors = {}
for name, method in [("Forward Euler", partial(ode_rk_explicit, table=euler)),
                     ("Heun", partial(ode_rk_explicit, table=heun)),
                     ("RK4", partial(ode_rk_explicit, table=rk4))]:
    errors = []
    for h in hs:
        error = mms(eqn, method, u_0, h=h, T=8)
        errors.append(error)
    method_errors[name] = numpy.asarray(errors)
    pyplot.loglog(hs, errors, ".", label=name)
    
pyplot.loglog(hs, 1.5*hs, label="$\mathcal{O}(h)$")
pyplot.loglog(hs, 2*hs**2, label="$\mathcal{O}(h^2)$")
pyplot.loglog(hs, 0.1*hs**4, label="$\mathcal{O}(h^4)$")
pyplot.legend()
pyplot.xlabel("h")
pyplot.ylabel("$e_{n,h}$");
```

### Questions

1. What timestep would you need for Heun's method and forward Euler respectively to obtain an error of $10^{-10}$?


```python

```

### Work precision diagrams

This comparison is slightly unfair, since when we are only looking at the error, higher-order methods will clearly perform better for the same timestep. Especially on this oscillatory test problem that has a very smooth solution.

A fairer method of comparison is to use _work precision_ diagrams. In these, we plot some proxy for the total _cost_ of a method against the obtained error. For explicit methods, the proxy is usually the _number of function evaluations_ (since it avoids arguments about inefficiencies in implementation) though _time to solution_ is also valid.

With this in mind, we can replot our convergence data.


```python
def work(hs, table, T=15):
    A = table.A
    b = table.b
    stages, _ = A.shape
    if numpy.allclose(A[-1, :], b):
        stages -= 1
    steps = T/hs
    return steps * stages
```


```python
cost = partial(work, hs, T=15)
pyplot.figure()
for name, table in [("Forward Euler", euler),
                    ("Heun", heun),
                    ("RK4", rk4)]:
    pyplot.loglog(cost(table), method_errors[name], marker="o", label=name)
pyplot.legend()
pyplot.xlabel("Number of function evaluations")
pyplot.ylabel("Error")
pyplot.title("Error vs Cost");
```

### Questions

1. This looks like the high-order method still always wins. But can you actually increase the step size sufficiently that you can extrapolate to the low-error case?


```python

```

## Error control

In CMIII you saw a relatively simple method for adaptive timestepping. Given a single timestepping scheme with step length $h$, the error can be estimated by performing two steps with step length $h/2$ and computing the difference in the resulting guess.

This is rather wasteful, in that we need to perform three time integrations to take one step.

For higher-order RK schemes, a more sophisticated approach is possible. Since the constraints on the coefficients in the Butcher table can admit more than one possible solution it is often possible to use the same set of stages with two different sets of *completion weights* (the $b$ vectors) to produce two updates of different order effectively "for free". The error can then be estimated by computing the difference between the two steps.

### Error estimation with embedded methods

Let us suppose we have an RK method

$$
\left[
\begin{array}{c|c}
c & A\\
\hline
& b^T\\
& \tilde{b}^T
\end{array}
\right]
$$

Where the completion weights $b^T$ produce a method of order $p$, and the weights $\tilde{b}^T$ produce a mthod of order $\tilde{p} = p - 1$. We then have, from a starting point $u_n$, two candidate solutions:

$$
\begin{align}
u_{n+1} &= u_n + h b^T f(Y)\\
\tilde{u}_{n+1} &= u_n + h \tilde{b}^T f(Y)\\
\end{align}
$$

where $f(Y)$ is the vector of intermediate function evaluations.

We can now estimate the error in a single step by

$$
e(h) = \| u_{n+1} - \tilde{u}_{n+1} \| = \| h (b - \tilde{b})^T f(Y) \| \in \mathcal{O}(h^p)
$$

Given some desired tolerance $\epsilon$, we want to find $h_*$ such that

$$
e(h_*) < \epsilon
$$

Since $e(h) \in \mathcal{O}(h^p)$ we can write $e(h) = c h^p$, then our requirement on $h_*$ means we want

$$
c h_*^p < \epsilon
$$

and so

$$
h_* < \left(\frac{\epsilon}{c}\right)^{1/p}.
$$

We estimate the unknown $c$ by using our error estimate:

$$
e(h) = c h^p \Leftrightarrow c = \frac{e(h)}{h^p}
$$

and substituting in

$$
h^* < h \left(\frac{\epsilon}{e(h)}\right)^{1/p}
$$

### BS3

A low-order method which implements this embedded idea is the [Bogacki-Shampine](https://en.wikipedia.org/wiki/Bogacki–Shampine_method) method of order 3, which comes with an embedded 2nd order method.


```python
bs3 = ButcherTable([[0, 0, 0, 0],
                    [1/2, 0, 0, 0],
                    [0, 3/4, 0, 0],
                    [2/9, 1/3, 4/9, 0]],
                   [2/9, 1/3, 4/9, 0], # order 3
                   order=3,
                   btilde=[7/24, 1/4, 1/3, 1/8]) # order 2
```

Although this method has 4 stages, it only requires 3 function evaluations, since if we stare closely we can see that the last stage has the same coefficients as the order 3 completion formula (so we can reuse the last stage for the next timestep). This is termed the _First Same as Last_ or FSAL property.


```python
bs3
```

### ODE45

The default ODE solver in MATLAB (amongst others) is the 5th order [Dormand-Prince](https://en.wikipedia.org/wiki/Dormand–Prince_method) method (with a 4th order embedded method). One can vaguely imagine that this table was not computed by frantically Taylor expanding!


```python
dp45 = ButcherTable([[0, 0, 0, 0, 0, 0, 0],
                     [1/5, 0, 0, 0, 0, 0, 0],
                     [3/40, 9/40, 0, 0, 0, 0, 0],
                     [44/45, -56/15, 32/9, 0, 0, 0, 0],
                     [19372/6561, -25360/2187, 64448/6561, -212/729, 0, 0, 0],
                     [9017/3168, -355/33, 46732/5247, 49/176, -5103/18656, 0, 0],
                     [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84, 0]],
                    [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84, 0], 
                    order=5,# 5th order
                    btilde=[5179/57600, 0, 7571/16695, 393/640, -92097/339200, 187/2100, 1/40]) # 4th order
dp45
```

Note how this method also has the FSAL property for the 5th order completion formula.

## Adaptive error control

Let us look at an example use of adaptive error control. The step-size acceptance here is rather simple, for a more thorough treatment, see [Hairer & Wanner (1993), pages 167-168](https://www.springer.com/gp/book/9783540566700).

We solve the coupled chemical reaction system known as the "Brusselator" [Lefever & Nicolis (1971)](https://doi.org/10.1016/0022-5193(71)90054-3).

$$
\begin{align}
\dot{u}_1 &= A + u_1^2 u_2 - (B+1)u_1\\
\dot{u}_2 &= Bu_1 - u_1^2 u_2\\
A &= 1\\
B &= 3\\
u_1(0) &= 1.5\\
u_2(0) &= 3.0\\
\end{align}
$$


```python
class brusselator(object):
    def __init__(self,  A, B):
        self.A = A
        self.B = B
    
    def __call__(self, t, u):
        u1, u2 = u
        A = self.A
        B = self.B
        return numpy.asarray([A + u1**2*u2 - (B + 1)*u1, B*u1 - u1**2*u2])
```

For the error control we will use the embedded method in our RK scheme to calculate an error $\|u - \tilde{u}\|$. If this is larger than some specified $\epsilon$ we should reject the step. Additionally, an optimal $h$ for the next step (irrespective of whether we accepted the step or not) is chosen by using the inequality above.

We set:

$$
h_\text{next} = 0.9 * h_\text{current} \left( \frac{\| u - \tilde{u}\|}{\epsilon}\right)^{-1/p}
$$

where $p$ is the order of the method. The factor of 0.9 is a safety factor to ensure we have some wiggle room in our optimal step size estimation.


```python
def ode_rk_explicit_adaptive(f, u_0, table, eps=0.1, h_0=0.1, T=1):
    """Integrate f using an explicit RK scheme and adaptive timestepping.
    
    :arg f: The RHS a callable f(t, u) -> udot
    :arg u_0: initial value for u
    :arg table: a :class:`ButcherTable`.
    :arg eps: target local error.
    :arg h_0: Initial step size.
    :arg T: Final time."""
    A = table.A
    b = table.b
    btilde = table.btilde
    if btilde is None:
        raise ValueError("Provided scheme %r does not have embedded scheme" % table)
    order = table.order
    c = A.sum(axis=1)
    N = len(c)
    t = 0
    u_0 = numpy.asarray(u_0)
    u = u_0.copy()
    hist = [(t, u.copy(), h_0)]
    fY = numpy.zeros(u_0.shape + (N, ))
    h = h_0
    while t < T:
        accept = True
        h = min(h, T - t)
        fY[:] = 0
        # Standard RK stages
        for i in range(N):
            Yi = u.copy()
            for j in range(i):
                Yi += h * A[i, j] * fY[:, j]
            fY[:, i] = f(t + h*c[i], Yi)  
        # Error estimation
        ulo_up = h * fY @ btilde
        u_up = h * fY @ b
        err = numpy.linalg.norm(u_up - ulo_up)
        relerr = err / eps
        if relerr > 1: # error bigger than eps
            # The step was too big
            accept = False
        if relerr > 0:
            next_h = h * 0.9*relerr**(-1/order)
        else:
            # Some large number
            next_h = h*2
        if accept:
            # Step was fine (error was controlled), so use it.
            u += u_up
            t += h
            hist.append((t, u.copy(), h))
        # New step
        h = next_h
    ts, us, hs = zip(*hist)
    return numpy.asarray(ts), numpy.asarray(us), numpy.asarray(hs)
```


```python
u = numpy.asarray([1.5, 3.0])
eqn = brusselator(A=1, B=4)

thist, uhist, hhist = ode_rk_explicit_adaptive(eqn, u, bs3, eps=1e-3, h_0=0.1, T=20)
```


```python
pyplot.figure()
pyplot.plot(thist, uhist);
pyplot.figure()
pyplot.plot(thist, hhist);
pyplot.xlabel("t")
pyplot.ylabel("h");
```

## Questions

1. Use a small timestep, high-order integrator to compute an "exact" solution. Use this to produce MMS convergence plots for the BS3 and DP45 methods using this adaptive error control mechanism. How do you need to vary $\epsilon$ to get the convergence order you expect?

2. Using the same high-order "exact" solution, produce a work-precision diagram comparing fixed step-size methods against adaptive methods. Are the adaptive methods always superior in terms of work for a given accuracy? You may find it useful to augment the `brusselator` class with a count object that keeps track of the number of function evaluations (don't forget to create a fresh object for each integration though!).

3. Compare, using a work-precision diagram, the cost of this "embedded" adaptivity scheme with the "try with h/2 and compare" approach of CMIII.


```python

```

# Back to stiff problems

## Implicit Runge-Kutta methods

We have already seen that we will want an implicit, rather than explicit, method when taking large timesteps. The same formalism for expressing Runge-Kutta methods is also applicable to implicit ones. In this case, $A$ is no longer strictly lower-triangular. Popular implicit RK schemes are typically at least lower triangular ($a_{i,j} = 0, j > i$). This way, although we need to solve a system for each stage, we can still do forward-substitutions, and the stage equations become:

$$
Y_i - h a_{i,i} f(t_n + c_i h, Y_i) = u_n + h \sum_{j < i} a_{i, j} f(t_n + c_j h, Y_j).
$$


These are termed DIRK (_diagonally_ implicit Runge-Kutta) methods. 

A subclass that is even better (if achievable) is SDIRK (_singly_ diagonally implicit Runge-Kutta) methods. This last implicit class has $a_{i,i} = \alpha$ for all $i$.

These methods are attractive because we can often reuse (part of) the setup cost for solving the systems at each stage.

Let's return to the _stiff_ equation we first saw

$$
\dot{u} = -k(u - \cos t)
$$

and try solving it when $k$ is very large (this increases the stiffness).



```python
class stiff(object):
    def __init__(self, k):
        self.k = k
    def f(self, t, u):
        return -self.k * (u - numpy.cos(t))
    
    def explicit(self, t, u):
        # This part of the forcing function will be moved to the rhs
        return self.k * numpy.cos(t)
    
    __call__ = f
    
    def J(self, t, u):
        # This is the linearly implicit part of the forcing function (solved for)
        # df/du
        return numpy.asarray(-self.k).reshape(1, 1)
    
    def exact(self, t, u_0):
        k = self.k
        k2p1 = k/(1 + k**2)
        return (u_0 - k*k2p1)*numpy.exp(-k*t) + k2p1 * (numpy.sin(t) + k*numpy.cos(t))
```


```python
def ode_dirk_linear(eqn, u_0, table, h=0.1, T=1):
    A = table.A
    b = table.b
    order = table.order
    c = A.sum(axis=1)
    N = len(c)
    for i in range(N):
        if any(A[i, i+1:] > 0):
            raise ValueError("Only for diagonally implicit RK schemes")
    t = 0
    u_0 = numpy.asarray(u_0).reshape(-1)
    u = u_0.copy()
    hist = [(t, u.copy())]
    fY = numpy.zeros(u_0.shape + (N, ))
    while t < T:
        h = min(h, T - t)
        fY[:] = 0
        # Standard RK stages
        for i in range(N):
            Yi = u.copy()
            tstage = t + h*c[i]
            for j in range(i):
                Yi += h * A[i, j] * fY[:, j]
            if A[i, i] != 0:
                # Solving
                # Yi - h a_ii f(tstage, Yi) = u_n + \sum_{i<j} h A[i, j] f(t + c[j] h, Y_j)
                # f(stage, Yi) splits into linearly implicit part and (non)linear explicit part (moved to rhs)
                J = eqn.J(tstage, Yi)
                rhs = Yi + h * A[i, i] * eqn.explicit(tstage, Yi)
                Yi = numpy.linalg.solve(numpy.eye(J.shape[0]) - h * A[i, i] * J, rhs)
            fY[:, i] = eqn(tstage, Yi)
        u += h * fY @ b
        t += h
        hist.append((t, u.copy()))
    ts, us = zip(*hist)
    return numpy.asarray(ts), numpy.asarray(us)
```


```python
beuler = ButcherTable([[1]], [1], order=1)

test = stiff(200)

u_0 = numpy.array(0.2)
thist, uhist = ode_dirk_linear(test, u_0, beuler, h=0.1, T=12)
pyplot.figure()
pyplot.plot(thist, uhist, "o", linestyle="solid", label='Backward Euler')
pyplot.plot(thist, test.exact(thist, u_0), label='exact');
pyplot.legend();
```


```python
# 4 stage, 3rd order scheme from Ascher et al. (1997) https://doi.org/10.1137/0732037
sdirk43 = ButcherTable([[1/2, 0, 0, 0],
                       [1/6, 1/2, 0, 0],
                       [-1/2, 1/2, 1/2, 0],
                       [3/2, -3/2, 1/2, 1/2]],
                      [3/2, -3/2, 1/2, 1/2], order=3)
sdirk43
```


```python
u_0 = numpy.array(0.2)

test = stiff(k=1000)
thist, uhist = ode_dirk_linear(test, u_0, sdirk43, h=0.1, T=20)
pyplot.figure()
pyplot.plot(thist, uhist, "o", linestyle="solid", label='DIRK43')
pyplot.plot(thist, test.exact(thist, u_0), label='exact');
pyplot.legend()
```

## Questions

1. Produce work-precision diagrams for the implicit methods on the stiff equation above. Compare them to using an explicit integrator instead.


```python

```

2164fcd680e221414b1edfeb962b2a9aba52798e
ipynb
Jupyter Notebook
material/time-integration-stability-rk.ipynb
wence-/comp4031-scicomp
2840ebb8e1cb09b865400c0639b247005fa3c601
BSD-2-Clause
2019-10-07T10:43:06.000Z
2019-10-07T10:43:06.000Z
material/time-integration-stability-rk.ipynb
wence-/comp4031-scicomp
2840ebb8e1cb09b865400c0639b247005fa3c601
BSD-2-Clause
material/time-integration-stability-rk.ipynb
wence-/comp4031-scicomp
2840ebb8e1cb09b865400c0639b247005fa3c601
BSD-2-Clause
2019-11-28T13:26:48.000Z
2019-11-28T13:26:48.000Z
Qwen/Qwen-72B
1. YES
2. YES