# Spatial discretisation

So far, we've seen time derivatives and ordinary differential equations of the form

$$
\dot{u} = f(t, u).
$$

Most problems one encounters in the real world have spatial as well as time derivatives. Our first example is the [*Poisson equation*](https://en.wikipedia.org/wiki/Poisson%27s_equation):

$$
\begin{align}
-\frac{\text{d}^2 u}{\text{d} x^2} &= f(x) \quad x \in \Omega = (-1, 1)\\
 u(-1) &= a\\
 \frac{\text{d} u}{\text{d} x}(1) &= b\\
 \end{align}.
$$

This is termed a *boundary value problem* (BVP), as opposed to the ODEs which were *initial value problems*, since we do not specify an initial condition, but rather a condition on the boundary of the domain.

This equation appears in a remarkably large number of places. As one example, it models the equilibrium temperature profile of a thermally conducting material maintained at constant temperature $a$ on the left and cooled at constant rate $b$ on the right.

To solve this problem numerically we must make a number of choices:

- how to represent the solution $u$;
- how to compute its derivatives;
- how to enforce the boundary conditions;
- how and where to evaluate $f$;
- in what sense we would like our solution to satisfy the equation.

[Iserles' book](http://www.damtp.cam.ac.uk/user/ai/Arieh_Iserles/Textbook.html) contains a quite mathematical treatment of finite differences. I also like [Randy LeVeque's *Finite difference methods for ordinary and partial differential equations*](http://staff.washington.edu/rjl/fdmbook/).

## Finite difference framework

We will focus on _finite difference_ methods here, which make the following choices in answer to the questions above:

- The solution $u(x)$ is represented by _pointwise_ values $u_i = u(x_i)$ at some discrete set of points $-1 = x_0 < x_1 < \dots < x_N = 1$. Importantly, the framework _does not_ specify the value of $u$ outside of these points;
- derivatives of $u$ at points $x_i$ are approximated using differencing formulae that utilise a finite number of neighbouring points (independent of $N$);
- boundary conditions are either enforced pointwise (e.g. the $u(-1)$ case above), or (when constraining derivatives) with one-sided differencing formulae;
- $f$ is evaluated pointwise at each $x_i$;
- we require that our finite difference method satisfies the equation pointwise at each $x_i$ in the interior of the domain.

### Differencing formulae

Our starting point is the definition of a derivative:

$$
\frac{d u(x)}{d x} = \lim_{\epsilon \to 0} \frac{u(x + \epsilon) - u(x)}{\epsilon}.
$$

If we wish to approximate this in our finite difference framework, where we only have point values, we can do so using neighbouring values. Writing $x_{i+1} - x_i = h$ for convenience, we can write

$$
\frac{d u(x_i)}{d x} \approx \frac{u(x_{i+1}) - u(x_i)}{h} = \frac{u_{i+1} - u_i}{h} =: D_+ u_i.
$$

This is a _one-sided_ approximation: we only use $u_i$ and $u_{i+1}$. Another one-sided approximation would be to offset in the other direction

$$
D_{-} u_i := \frac{u_i - u_{i-1}}{h}.
$$

Finally, we can also use a _centred_ approximation by averaging the two one-sided approximations:

$$
D_0 u_i := \frac{u_{i+1} - u_{i-1}}{2h} = \frac{1}{2} (D_+ + D_-) u_i.
$$

Let's have a picture


```python
%matplotlib notebook

import numpy
from matplotlib import pyplot
import matplotlib.lines as mlines
pyplot.style.use('ggplot')

n = 200
h = 2/(n-1)
x = numpy.linspace(1,2.5,n)
pyplot.plot(x, numpy.sin(x));

def newline(p1, p2, **kwargs):
    ax = pyplot.gca()
    xmin, xmax = ax.get_xbound()

    if(p2[0] == p1[0]):
        xmin = xmax = p1[0]
        ymin, ymax = ax.get_ybound()
    else:
        ymax = p1[1]+(p2[1]-p1[1])/(p2[0]-p1[0])*(xmax-p1[0])
        ymin = p1[1]+(p2[1]-p1[1])/(p2[0]-p1[0])*(xmin-p1[0])

    l = mlines.Line2D([xmin,xmax], [ymin,ymax], **kwargs)
    ax.add_line(l)
    return l

h = 0.25
xi = 1.6
ximinus = xi - h
xiplus = xi + h

pyplot.plot([ximinus, xi, xiplus], numpy.sin([ximinus, xi, xiplus]), marker="o", linestyle="none")

newline((xi, numpy.sin(xi)), (xiplus, numpy.sin(xiplus)), linestyle="dashed", label="$D_+ u(x)$")
newline((xi, numpy.sin(xi)), (ximinus, numpy.sin(ximinus)), linestyle="dotted", label="$D_- u(x)$")
newline((ximinus, numpy.sin(ximinus)), (xiplus, numpy.sin(xiplus)), linestyle="-.", label="$D_0 u(x)$")


newline((xi, numpy.sin(xi)), (xiplus, numpy.sin(xi) + h*numpy.cos(xi)), color="black", label="$u'(x)$")

pyplot.legend();
```

Let's look at the full approximate derivative too.


```python
n = 10
h = 2/(n-1)
x = numpy.linspace(-1,1,n)
u = numpy.sin(x)
pyplot.figure()
pyplot.plot(x, numpy.cos(x), label="$u'$");
pyplot.plot(x[1:], (u[1:] - u[:-1])/h, label="$D_-$", marker="o", linestyle="none")
pyplot.legend();
```

## Accuracy

Certainly from the picture of the slope above, it appears that the centered difference formula is more accurate than the one-sided approximations. Can we formalise this at all?

To do so, we turn to the favourite tool of the budding numericist, the *taylor expansion*.

### Recap, Taylor expansions

For a sufficiently smooth function $u$, given a point $x$, we can represent the function at a new point $x + h$ by its Taylor expansion

$$
u(x + h) = u(x) + u'(x) h + \frac{1}{2} u''(x) h^2 + \frac{1}{6} u'''(x) h^3 + \dots = \sum_{n=0}^\infty \frac{1}{n!} h^n u^{(n)}(x) 
$$

where the notation $u'$ is shorthand for $\frac{\text{d} u}{\text{d} x}$ and $u^{(n)}(x) = \frac{\text{d}^n u}{\text{d} x^n}$.

Let's look at how this works for a sample function.


```python
from functools import reduce
from operator import mul

def u(x, n=0):
    factor = (-1)**(n // 2)
    if n % 2 == 0:
        u_ = numpy.sin
    else:
        u_ = numpy.cos
    return factor * u_(x)

pyplot.figure()

x = numpy.linspace(0.5, 1.75, 500)

pyplot.plot(x, u(x));

x0 = 0.6
h = 0.8

def fac(n):
    return reduce(mul, range(1, n+1), 1)

def taylor(u, x0, h, n):
    return u(x0) + sum(h**i/fac(i) * u(x0, i) for i in range(1, n))

xs = numpy.linspace(x0, x0+h, 20)
for n in range(1, 5):
    pyplot.plot(xs, [taylor(u, x0, x - x0, n) for x in xs], marker="o", label=r"$\tilde{u} + \mathcal{O}(h^%d)$" % n)

pyplot.legend();
```

If we chop off the series at some finite $n$ we write

$$
u(x + h) = u(x) + u'(x) h + \frac{1}{2} u''(x) h^2 + \mathcal{O}(h^3)
$$

with $h$ sufficiently small.

To determine the order of a method, we substitute the Taylor expansion into the differencing expression and calculate.

As an example, let us consider the one-sided differencing operator $D_+$. To simplify notation, we will choose $x = 0$, and we have

$$
\begin{align}
u'(0) &\approx \frac{u(h) - u(0)}{h} \quad \text{ definition of } D_+\\
      &= h^{-1}(\underbrace{u(0) + u'(0) h + \frac{1}{2} u''(0) h^2 + \mathcal{O}(h^3)}_{u(h)} - u(0)) \\
      &= u'(0) + \frac{1}{2} u''(0) h + \mathcal{O}(h^2)\\
\end{align}.
$$

The leading-order error term in the right hand side is $\mathcal{O}(h)$, and so we say that this is a _first-order_ method. Derivation that the operator $D_-$ is also first-order proceeds identically.

#### Questions

1. Show that the centered difference operator $D_0$ computes a second-order accurate derivative.


```python

```

## Stability

We will postpone mathematical discussion of stability for a while, and give an intuition for some potential problems. Let us first check that our implementation of differencing operators provides us with the expected (mathematical) convergence orders for a smooth function.


```python
def dplus(x, u):
    return x[:-1], (u[1:] - u[:-1])/(x[1:] - x[:-1])

def dminus(x, u):
    return x[1:], (u[1:] - u[:-1])/(x[1:] - x[:-1])

def center(x, u):
    return x[1:-1], (u[2:] - u[:-2])/(x[2:] - x[:-2])
```


```python
grids = 2**numpy.arange(3, 10)

def error(f, df, op):
    for n in grids:
        x = numpy.linspace(-1, 1, n)
        x, y = op(x, f(x))
        yield numpy.linalg.norm(y - df(x), numpy.inf)

pyplot.figure()
for op in [dplus, dminus, center]:
    pyplot.loglog(1/grids, list(error(numpy.sin, numpy.cos, op)), marker="o", linestyle="none", label=op.__name__)
    
pyplot.xlabel("Resolution ($h$)")
pyplot.ylabel("$l_\infty$ error in derivative")

pyplot.loglog(1/grids, 1/grids, label="$h$")
pyplot.loglog(1/grids, 1/grids**2, label="$h^2$")
pyplot.legend();
```

So both the one-sided differences are first-order accurate, whereas the centered difference is second-order accurate, as expected. One thing to be wary of, however, is using some of these approximations for functions that are "rough" on the grid scale.

We can make this question more precise by asking whether there are functions whose derivatives are non-zero, but for which our numerical approximations compute $u'(x_i) = 0$.

Let's try and contrive an example:


```python
x = numpy.linspace(-1, 1, 9)
xf = numpy.linspace(-1, 1, 100)

def f(x):
    return numpy.cos(1/2 + 4*numpy.pi*x)

def df(x):
    return -4*numpy.pi*numpy.sin(1/2 + 4*numpy.pi*x)
```


```python
pyplot.figure()
pyplot.plot(x, f(x), marker="o", label="coarse")
pyplot.plot(xf, f(xf), "-", label="fine")
pyplot.legend();
```

What about the derivatives?


```python
pyplot.figure()
for op in [dplus, dminus, center]:
    x_, y = op(x, f(x))
    pyplot.plot(x_, y, "o-", label=op.__name__)
pyplot.plot(xf, df(xf), "-", label="Exact")
pyplot.legend();
```

The centered difference approximation produces a _zero_ derivative for this function. Hence if we have a solution $u(x)$, we can (at least to the numerical operator) construct a new solution $\tilde{u}(x) = u(x) + f(x)$.

Suddenly, even if our actual equation has a unique solution, the numerical solution is no longer unique. This turns out to cause all kinds of terrible issues with numerical algorithms and must be avoided at all costs.

## Higher-order derivatives

We can compute high-order derivatives by repeatedly applying differencing operators for lower-order derivatives.

For example, the second derivative

$$
\frac{\text{d}^2 u}{\text{d} x^2} \approx D^2 u_i = D_+ D_- u_i = \frac{1}{h^2}\left(u_{i+1} - 2 u_i + u_{i-1}\right) = D_- D_+ u_i.
$$

### Questions

1. Show that this is a _second-order_ accurate approximation of the second derivative
2. We could also use $D^2 u_i = D_0 D_0 u_i$, derive the stencil for this case.
3. Finally, show that if we define a "half-step" centered difference operator

$$
\hat{D}_0 u = \frac{1}{h}\left[u\left(x + \frac{h}{2}\right) - u\left(x - \frac{h}{2}\right)\right]
$$
then we have

$$
D^2 = \hat{D}_0 \circ \hat{D}_0 = D_+ \circ D_-
$$


```python

```

## Boundary conditions

The final missing piece required before we can solve the our first PDE is to figure out how we will treat boundary conditions. To do this, we will first recast the differencing operators as matrices. it is then somewhat easier to see what is going on.

We can think of the differencing operator acting on an entire vector

$$
U = \begin{bmatrix}
u_0\\
u_1\\
\vdots\\
u_N
\end{bmatrix}
$$

at once. For example, we can write

$$
D_+ = \frac{1}{h} \begin{bmatrix}
-1 & 1 & 0 & \dots & 0\\
0 & -1 & 1 & \dots & 0\\
\vdots & \ddots & \ddots & \vdots & \vdots\\
0 & \dots & 0 & -1 & 1\\
0 & \dots & 0 & 0 & -1\\
\end{bmatrix}
$$

and

$$
D^2 = \frac{1}{h^2} \begin{bmatrix}
-2 & 1 & 0 &  \dots & 0\\
1 & -2 & 1 &  \dots & 0\\
\vdots & \ddots & \ddots & \vdots & \vdots\\
0 & \dots & \dots 1 & -2 & 1\\
0 & \dots & 0 & 1 & -2\\
\end{bmatrix}.
$$

Recall, our problem was to find $u \in (-1, 1)$ satisfying

$$
\begin{align}
\frac{\text{d}^2 u}{\text{d} x^2} u &= f\\
u(-1) &= a\\
\frac{\text{d} u}{\text{d} x}(1) &= b
\end{align}
$$

In matrix form, this becomes

$$
\underbrace{D^2}_{A} \underbrace{\begin{bmatrix}
u_0\\
u_1\\
\vdots\\
u_N
\end{bmatrix}}_{U} = \underbrace{\begin{bmatrix} f_0\\f_1\\\vdots\\f_N\end{bmatrix}}_{F}.
$$

This works perfectly in the interior of the domain, but we need to figure out what to do at the boundaries. For example, we can't use the standard differencing operator for $D^2$ on $u_N$, because $u_{N+1}$ does not exist.
Fortunately, our boundary conditions inform how to modify the matrix appropriately.

### Dirichlet conditions

These conditions, of the form 

$$
u(-1) = a
$$

specify the _value_ of the solution at a particular point (or set of points). This means, that rather than solving a small equation to determine the value at this point, we _already know_ and can instead replace the relevant rows of the matrix. Let us suppose we have ordered our points such that $u(-1)$ corresponds to $u_0$. Then we have

$$
\begin{bmatrix}
1 & 0 & \dots & 0\\
 & & & \\
 & & A_{1:,:} & \\
 & & &\\
\end{bmatrix}
U = \begin{bmatrix} a\\ \\ F_{1:} \\ \\ \end{bmatrix}.
$$

In general, if we have a boundary value $\alpha_i$ that constrains $u_i$ then we replace the $i$th row with the identity, and the $i$th value in the right hand side with $\alpha_i$.

This modification destroys any symmetry that might have existed in the matrix $A$, since we have zeroed rows, but not the corresponding columns. If we write the linear system in block form, we can, however, see a way around this:

$$
\begin{bmatrix}
I & 0\\
A_{10} & A_{11}
\end{bmatrix}
\begin{bmatrix}
U_0\\
U_1
\end{bmatrix}
=
\begin{bmatrix}
F_0\\
F_1
\end{bmatrix}
$$

Since we know $U_0$ (they are just $F_0$), we can forward-substitute and move the lower-left block of the matrix onto the right hand side, to produce

$$
\begin{bmatrix}
I & 0\\
0 & A_{11}
\end{bmatrix}
\begin{bmatrix}
U_0\\
U_1
\end{bmatrix}
=
\begin{bmatrix}
F_0\\
F_1 - A_{10}F_0
\end{bmatrix}.
$$

This is often a convenient form to work with.

Alternately, since the equations for $U_0$ are just the identity, we can write our solver to just handle

$$
A_{11} U_1 = F_1 - A_{10}F_0
$$

and insert the boundary values into a big vector whenever we need to visualise it.

## Neumann conditions

We can now treat boundary conditions that constrain the value of the solution, but recall that the condition at $x=1$ instead constrains gradient of the solution. We cannot do this by setting values, but must instead form an equation for the boundary value.

### One-sided difference

There are typically two ways to do this. We either come up with a one-sided differencing formula for the derivative directly. For example, recalling the one-sided difference we might replace the boundary term

$$
\frac{\text{d} u}{\text{d} x}(1) = b
$$

by

$$
\frac{u_n - u_{n-1}}{h} = b.
$$

This is simple but has some potential drawbacks

1. We need to make a different choice for the discretisation on the boundary to that in the interior
2. This choice may not have the same order of accuracy
3. It may destroy symmetry that previously existed in the problem.

### Ghost values

An alternate option is to introduce a (or possibly more than one) _ghost value_ outside of the domain such that we can then just use our interior discretisation. We then define the value of this ghost point to be the reflection (possibly weighted by the boundary value) of the interior point.

That is, we introduce $u_{n+1} = u(x_{n+1})$ and set

$$
u_{n+1} = u_{n-1} + 2b(\underbrace{x_n - x_{n-1}}_{h}).
$$

Now we can use our interior discretisation

$$
\frac{-u_{n-1} + 2u_n - u_{n+1}}{h^2} = f(x_n)
$$

substituting in the definition of $u_{n+1}$ we obtain

$$
\begin{align}
\frac{-u_{n-1} + 2u_n - (u_{n-1} + 2bh)}{h^2} &= f(x_n)\\
\frac{2(u_n - u_{n-1})}{h^2} &= f(x_n) + \frac{2b}{h}\\
\frac{u_n - u_{n-1}}{h^2} &= \frac{f(x_n)}{2} + \frac{b}{h}
\end{align}.
$$

Let's compare these approaches.


```python
def laplacian(N, rhsfunc):
    x = numpy.linspace(0, 1, N+1)
    h = 1/N
    rhs = rhsfunc(x)
    e = numpy.ones(N)
    # interior discretisation
    L = (2*numpy.eye(N+1) - numpy.diag(e, 1) - numpy.diag(e, -1)) / h**2
    return x, L, rhs, h
    
def apply_dirichlet(L, rhs, h, vals, indices):
    N, _ = L.shape
    diag = numpy.eye(1, N)
    bcmask = numpy.zeros(N, dtype=bool)
    bcmask[indices] = True
    # Dirichlet rows
    L[numpy.ix_(bcmask)] = numpy.vstack([numpy.roll(diag, i) for i in indices])
    rhs[numpy.ix_(bcmask)] = vals
    # Forward substitute
    rhs[numpy.ix_(~bcmask)] -= L[numpy.ix_(~bcmask, bcmask)] @ vals
    L[numpy.ix_(~bcmask, bcmask)] = 0
    return L, rhs, h

def apply_neumann_oneside(L, rhs, h, b, index):
    N, _ = L.shape
    assert index == N - 1
    L[index, :] = 0
    L[index, index] = 1/h
    L[index, index - 1] = -1/h
    rhs[index] = b
    return L, rhs, h

def apply_neumann_ghost(L, rhs, h, b, index):
    N, _ = L.shape
    L[index, index] /= 2
    rhs[index] = b/h + rhs[index]/2
    return L, rhs, h
```

We'll solve

$$
\begin{align}
-\frac{\text{d}^2 u}{\text{d} x^2} &= e^x \text{ in } (0, 1)\\
u(0) &= e^{0}\\
\frac{\text{d} u}{\text{d} x} &= e^{1}\\
\end{align}
$$

with convenient exact solution $u(x) = e^{x}$.


```python
N = 10
rhsfunc = lambda x: -numpy.exp(x)
exact = lambda x: numpy.exp(x)
x, L, rhs, h = laplacian(N, rhsfunc)
L, rhs, h = apply_dirichlet(L, rhs, h, [exact(0)], [0])
L, rhs, h = apply_neumann_oneside(L, rhs, h, exact(1), N)
uoneside = numpy.linalg.solve(L, rhs)

x, L, rhs, h = laplacian(N, rhsfunc)
L, rhs, h = apply_dirichlet(L, rhs, h, [exact(0)], [0])
L, rhs, h = apply_neumann_ghost(L, rhs, h, exact(1), N)
ughost = numpy.linalg.solve(L, rhs)
```


```python
pyplot.figure()
pyplot.plot(x, uoneside, label="Computed Oneside")
pyplot.plot(x, ughost, label="Computed Ghost")
pyplot.plot(x, exact(x), label="Exact")
pyplot.legend();
```

## Observations

Perhaps unsurprisingly, the one-sided application of the Neumann conditions performs worse than the ghost version. Interestingly, they are effectively the *same* discretisation in the matrix, the only difference is that in the ghost version, we corrected the right hand side we're solving for by a small amount to take into account the issues in the one-sided discretisation.

It looks like we have a lower-order scheme. Let's check by performing an MMS test.

To do so, we have to introduce how to measure errors. Since our discrete solution $u_i$ is supposed to approximate $u(x_i)$ it is natural to consider the pointwise errors $u_i - u(x_i)$. Let us now consider how to measure the size of the error vector (or indeed any vector).

$$
E = \begin{bmatrix} u_0 - u(x_0)\\
\vdots\\
u_n - u(x_n)
\end{bmatrix}
$$

So far we have been using the $\infty$-norm or $\max$-norm.

$$
\|E\|_\infty := \max_{0 \le i \le n} |E_i| = \max_{0 \le i \le n} |u_i - u(x_i)|
$$

which measures the largest pointwise error over the interval.

Other common norms are the $1$-norm

$$
\|E\|_1 = h \sum_{i=0}^n |E_i|
$$

and the $2$-norm

$$
\|E\|_2 = \left(h \sum_{i=0}^n |E_i|^2 \right)^{1/2}.
$$

Notice the factor of $h$ appearing in these definitions. This is needed so the norm does not spuriously grow when we add more points.

### Aside

These are special cases of $l_p$ norms

$$
\|E\|_p = \left(h \sum_{i=0}^n |E_i|^p\right)^{1/p}.
$$


```python
# Here we use the 2-norm
def error(u, exact, h):
    return numpy.sqrt(h)*numpy.linalg.norm(u - exact)
```


```python
def mms_errors(neumann):
    errors = []
    Ns = numpy.asarray(list(2**i for i in range(4, 11)))
    rhsfunc = lambda x: -numpy.exp(x)
    exact = lambda x: numpy.exp(x)
    for N in Ns:
        x, L, rhs, h = laplacian(N, rhsfunc)
        L, rhs, h = apply_dirichlet(L, rhs, h, [exact(0)], [0])
        L, rhs, h = neumann(L, rhs, h, exact(1), N)
        u = numpy.linalg.solve(L, rhs)
        errors.append(error(u, exact(x), 1/N))
    return 1/Ns, numpy.asarray(errors)
```


```python
_, oneside = mms_errors(apply_neumann_oneside)
hs, ghost = mms_errors(apply_neumann_ghost)
pyplot.figure()
pyplot.loglog(hs, oneside, "o", label="Oneside");
pyplot.loglog(hs, ghost, "x", label="Ghost");
pyplot.loglog(hs, hs, label="$\mathcal{O}(h)$");
pyplot.loglog(hs, hs**2, label="$\mathcal{O}(h^2)$");
pyplot.legend();
```

This confirms our suspicion that the one-sided differencing for the Neumann condition is only first-order accurate.

An alternative approach to obtaining a second-order scheme (rather than the ghost method above) is to try and determine a second-order accurate one-sided difference approximation to the first derivative. We will state an example first, and then see where it comes from.

A second-order accurate one-sided approximation to the first derivative is obtained with

$$
\frac{\text{d} u}{\text{d} x} \approx \frac{1}{h}\left(\frac{3}{2} u_i - 2 u_{i-1} + \frac{1}{2} u_{i-2}\right)
$$


```python
def apply_neumann_oneside_second(L, rhs, h, b, index):
    N, _ = L.shape
    assert index == N - 1
    L[index, :] = 0
    L[index, index] = 3/(2*h)
    L[index, index - 1] = -2/h
    L[index, index - 2] = 1/(2*h)
    rhs[index] = b
    return L, rhs, h
```


```python
_, oneside = mms_errors(apply_neumann_oneside)
_, second = mms_errors(apply_neumann_oneside_second)
hs, ghost = mms_errors(apply_neumann_ghost)
pyplot.figure()
pyplot.loglog(hs, oneside, "o", label="Oneside");
pyplot.loglog(hs, second, "s", label="Oneside second order");
pyplot.loglog(hs, ghost, "x", label="Ghost");
pyplot.loglog(hs, hs, label="$\mathcal{O}(h)$");
pyplot.loglog(hs, hs**2, label="$\mathcal{O}(h^2)$");
pyplot.legend();
```

This converges at second order as, perhaps, expected. The absolute error is a little worse than the ghosted version. However, this approach is sometimes more convenient, especially on irregularly spaced meshes.

## Deriving high-order finite difference stencils

Where did the approximation

$$
\frac{\text{d} u}{\text{d} x} \approx \frac{1}{h}\left(\frac{3}{2} u_i - 2 u_{i-1} + \frac{1}{2} u_{i-2}\right)
$$

come from? Given some points at which we're allowed to evaluate $u$, we can derive an appropriate formula from the Taylor series using the *method of undetermined coefficients*. This works in a very similar way to determining the truncation error for a given expansion.

For the example above, we want to approximation $u'(x)$ and we are given $u_i = u(x)$, $u_{i-1} = u(x - h)$, and $u_{i-2} = u(x - 2h)$. We can write our differencing operator as a *linear combination* of the provided points

$$
D_2 u(x) = a u(x) + bu(x-h) + c u(x - 2h)
$$

where our goal is to determine $a$, $b$, and $c$ to minimise the truncation error (that is, give the best possible accuracy).

Let's Taylor-expand on the right hand side

$$
D_2 u(x) = a u(x) + b \overbrace{\left(u(x) - hu'(x) + \frac{h^2}{2} u''(x) - \frac{h^3}{6} u'''(x)\right)}^{u(x - h)} + c \overbrace{\left(u(x) - 2h u'(x) + \frac{4h^2}{2}u''(x) - \frac{8 h^3}{6} u'''(x)\right)}^{u(x-2h)} + \mathcal{O}(h^4)
$$

gathering terms we have

$$
D_2 u(x) = (a + b + c)u(x) - (b + 2c) h u'(x) + \frac{1}{2}(b + 4c)h^2 u''(x) - \frac{1}{6}(b + 8c) h^3 u'''(x) + \mathcal{O}(h^4).
$$

To maximise the accuracy of agreement with $u'(x)$ we need

$$
\begin{aligned}
a + b + c &= 0 && \text{zeroing the $h^0 u(x)$ term}\\
b + 2c &= -\frac{1}{h} && \text{ensuring that we have a $u'(x)$ term}\\
b + 4c &= 0 && \text{zeroing the $h^2 u''(x)$ term}\\ 
b + 8c &= 0 && \text{zeroing the $h^3 u'''(x)$ term}.
\end{aligned}
$$

Since we have only three unknowns, we can only satisfy three equations. To maxmimise the accuracy, we'll choose to zero the $h^2$ and $h^0$ terms, and live with the $h^3$ term. We therefore need to solve the linear system

$$
\begin{bmatrix}
1 & 1 & 1\\
0 & 1 & 2\\
0 & 1 & 4
\end{bmatrix}
\begin{bmatrix}
a \\ b \\ c
\end{bmatrix}
= 
\begin{bmatrix}
0 \\ -\frac{1}{h} \\ 0
\end{bmatrix}.
$$


```python
import numpy
A = numpy.asarray([[1, 1, 1],
                   [0, 1, 2],
                   [0, 1, 4]])
b = numpy.asarray([0, -1, 0])

numpy.linalg.solve(A, b)
```

So we have

$$
\begin{bmatrix}
a \\ b \\ c
\end{bmatrix} =
\frac{1}{2h}
\begin{bmatrix}
3\\
-4\\
1
\end{bmatrix}
$$

and hence our optimal formula is

$$
D_2 u(x) = \frac{1}{2h}(3 u_i - 4 u_{i-1} + u_{i-2})
$$

as advertised. We can immediately determine the accuracy of this approximation since we know the first term we did not manage to match exactly is

$$
- \frac{1}{6}(b + 8c) h^3 u'''(x)
$$

substituting in the values for $b$ and $c$ we have

$$
\begin{aligned}
D_2 u(x) - u'(x) &= -\frac{1}{6}\left(\frac{-2}{h} + \frac{8}{2h}\right) h^3 u'''(x) + \mathcal{O}(h^4)\\
                 &= -\frac{1}{3} h^2 u'''(x) + \mathcal{O}(h^4)
\end{aligned}
$$

so this approximation is second order accurate.


```python
def dtwo(x, u):
    h = x[2:] - x[1:-1]
    du = 1/(2*h) * (3 * u[2:] - 4*u[1:-1] + u[:-2])
    return x[2:], du

grids = 2**numpy.arange(3, 10)

def error(f, df, op):
    for n in grids:
        x = numpy.linspace(-1, 1, n)
        x, y = op(x, f(x))
        yield numpy.sqrt(1/n)*numpy.linalg.norm(y - df(x), None)

pyplot.figure()
pyplot.loglog(1/grids, list(error(numpy.sin, numpy.cos, dtwo)), marker="o", linestyle="none", label=dtwo.__name__)
    
pyplot.xlabel("Resolution ($h$)")
pyplot.ylabel("$l_2$ error in derivative")

pyplot.loglog(1/grids, 1/grids**2, label="$h^2$")
pyplot.loglog(1/grids, 1/grids**2, label="$h^2$")
pyplot.legend();
```


```python
numpy.linalg.norm?
```

## Questions

1. Use this technique to derive a third-order accurate one-sided differencing operator for $\frac{\text{d}}{\text{d} x}$ using 4 points $u(x)$, $u(x - h)$, $u(x - 2h)$, $u(x - 3h)$.


```python

```

## Differencing for advection

For centered difference approximations, there is no asymmetry in the stencil. For one-sided approximations, however, there is. If we just consider a three-point region centered around a point $i$, we can write the three approximations to $\frac{\text{d}}{\text{d} x}$ we have seen as stencils:

$$
\begin{aligned}
D_+ &= \frac{1}{h}\begin{bmatrix}0 & -1 &1\end{bmatrix}\\
D_- &= \frac{1}{h}\begin{bmatrix}-1 &  1 & 0\end{bmatrix}\\
D_0 &= \frac{1}{2h}\begin{bmatrix}-1 & 0 & 1\end{bmatrix}\\
\end{aligned}.
$$

Recall that we noticed that the centered difference approximation sometimes gave catastrophic results (all zero derivatives) for very rough functions. The question therefore might arise how to pick between $D_+$ and $D_-$. We will study this using the linear *advection equation* as a prototype.

This equation models the transport of some material by a bulk motion. This (especially when talking about fluid flow) is called convection. As usual [wikipedia has lots of information](https://en.wikipedia.org/wiki/Advection).

It looks remarkably benign, find $u(x, t)$ satisfying

$$
\partial_t u + c \cdot \nabla u = f(t, x)
$$

Where $c$ is the advecting velocity and

$$
\nabla u = \partial_x u
$$
in one dimension, and
$$
\nabla u = \begin{bmatrix} \partial_x u\\ \partial_y u\end{bmatrix}
$$
in two dimensions.

This is a first-order PDE, for which we need to supply one boundary condition (to pin down the spatial derivative) and one initial condition (to start everything off). The boundary condition, it turns out, has to be at the *inflow* boundary.

Let's try and solve this equation with an explicit Euler time integration scheme on the interval $[0, 5]$ and look at the effect of the different differencing operators. If we write this out we have (using superscripts for time points and subscripts for spatial points)

$$
u^{n+1}_i = u^n_i + \Delta t (f^n_i - c D u^n_i).
$$

We'll set the forcing function $f$ to be zero and pick boundary condition

$$
u(t, 0) = 0
$$

and initial condition

$$
u(0, x) = e^{-2(x - 2.5)^2}
$$


```python
L = 30
nx = 50
x = numpy.linspace(0, L, nx + 1)
h = L/nx
u = numpy.exp(-2*(x - L/2)**2)
uhat = numpy.zeros_like(x)
uhat[nx//3:2*nx//3] = 1


def Aupwind(nx, h, c):
    A = numpy.zeros((nx+1, nx+1), dtype=float)
    for i in range(1, nx+1):
        A[i, i-1] = c/h
        A[i, i] = -c/h
    # Boundary condition (will be fixed later)
    A[0, 0] = 0
    return A

def Adownwind(nx, h, c):
    A = numpy.zeros((nx+1, nx+1), dtype=float)
    for i in range(1, nx):
        A[i, i] = c/h
        A[i, i+1] = -c/h
    # Boundary condition (will be fixed later)
    A[0, 0] = 0
    A[nx, nx] = c/h
    return A

t = 0
tfinal = 10
c = 1
dt = 2/(nx + 1)

Id = numpy.eye(nx+1)

downwind = Id + dt*Adownwind(nx, h, c)
upwind = Id + dt*Aupwind(nx, h, c)

u = uhat
hist = [(t, u)]
while t < tfinal:
    u = upwind @ u
    hist.append((t, u))
    t += dt
hist = numpy.asarray(hist)    
```


```python
pyplot.figure()
for t, hist_ in hist[::len(hist)//10]:
    pyplot.plot(x, hist_, label=f"$t = {t:3.1f}$")
pyplot.legend();
```

## Observations

1. The "downwind" discretisation is unstable (blowing up near the boundary). The upwind version is stable, if we have a small enough time step. This makes physical sense, because the downwind discretisation is trying to obtain information "from the future".

2. The upwind discretisation is *dissipative*: the correct physical solution is just to transport the initial condition to the right, but we see that the peak spreads and flattens.

3. This also occurs with sharp fronts (e.g. transporting a hat function)

4. Some of this gets easier if we naturally have some diffusion (a second derivative) in the equation (as in the coursework).


```python

```

194fb871d95d6f3af6a42baf68807d07798824f4
ipynb
Jupyter Notebook
material/finite-differences.ipynb
wence-/comp4031-scicomp
2840ebb8e1cb09b865400c0639b247005fa3c601
BSD-2-Clause
2019-10-07T10:43:06.000Z
2019-10-07T10:43:06.000Z
material/finite-differences.ipynb
wence-/comp4031-scicomp
2840ebb8e1cb09b865400c0639b247005fa3c601
BSD-2-Clause
material/finite-differences.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