Lecture 6 (a)

Semi-discretization (Method of Lines)

Unstructured grids

Most problems have complex domains.
Boundary conditions crucial.
Little “structure” to problem.
What use is “structure” in the grid?
Use unstructured grids.
But how, whilst keeping accuracy?

An unstructued grid representing parts of the UK

Semi-discrete Methods

Take a “general” PDE

\[ \partial_t \phi = −\partial_x f(\phi) + \partial_x ( \kappa \partial_x \phi) = \mathcal{F}(\phi). \]

Use eg finite differences to approximate \(\mathcal{F} \to F\):

\[ \partial_t \phi_j = F(\phi_j) \simeq \frac{f_{i+1/2} - f_{i-1/2}}{\Delta x} + \dots \]

These are coupled ODEs. Solve using (eg) Runge-Kutta methods.

One significant complexity with finite difference methods is the coupling of differencing in space and in time. Stability and accuracy often mean that different approaches are needed. It’s also clear that there’s very different geometric aspects in the two cases, with complex boundaries usually only an issue in space, not in time.

One approach that can simplify our life is semi-discretization. Here we introduce a grid in space alone. We can then approximate all spatial derivatives by combinations of values on the spatial grid. This is usually represented by a general operator notation, where all the spatial derivatives are encapsulated by the curly F operator, and the discrete approximation to it by the upright F operator. The curly F operator takes functions as input and returns a function as output. The discrete approximation to the operator takes the discrete values of the function at all points of the grid and returns an approximation to the appropriate combinations of derivatives as a number at a single point on the grid.

The key advantage of the semi-discrete approach is that we are now solving ODEs. We have decoupled the behaviour in space and time. This makes it easier to construct complex approximations, which is important for high order methods. The ODEs are themselves coupled (as the spatial derivatives and their approximations will couple neighbouring gridpoints) and nonlinear, but standard numerical approximations for solving ODEs (such as Runge-Kutta methods) are well developed and easy to employ.

The disadvantages of semi-discrete methods are manifold, but usually minor. Fully discrete approaches at the same order of accuracy are usually more efficient and accurate in the absolute sense. It is also often true that formal stability or accuracy analysis in space and time separately does not give the complete information about the properties of the full scheme. However, the ease with which higher order methods can be designed and implemented in complex geometries makes semi-discrete approaches extremely useful in practice.

Example: FTBS

\[ \begin{aligned} \partial_t \phi_j &= -u \left. \left( \partial_x \phi \right) \right|_{x_j} \\ & \to -u \left( \phi_j - \phi_{j-1} \right) / \Delta x \, . \end{aligned} \]

def dphidt_bs(phi, dx, u=1):
    dphidt = np.zeros_like(phi)
    dphidt[:-1] = -u * (phi[:-1] - phi[1:]) / dx
    dphidt[-1] = dphidt[1] # Periodic, 1 ghost point
    return dphidt
def euler_step(phi_n, dphidt, dx, dt, u=1):
    return phi_n + dt * dphidt(phi_n, dx, u)
...
t=0
while t < t_end:
    t += dt
    phi_n = euler_step(phi_n, dphidt_bs, dx, dt)

Here is a brief example of the steps.

First we rewrite the standard advection equation as “time derivative of unknowns equals F”, where \(F\) is the right hand side. We replace all spatial derivatives with a discrete approximation. Here that approximation is backward space differencing. This gives the set of ODEs.

The implementation in code then splits the update into two pieces. The first function implements the spatial differencing: as can be seen that only needs the spatial grid spacing dx. This returns the time derivative of every point on the grid, including the boundary conditions. The second function solves the ODE in time using the explicit Euler method, which is precisely what Forward Time differencing gives.

The final part of the code steps through the update within a loop to actually solve the problem.

This is likely many more lines of code than a simple, fully discrete implementation would be. However, we can see that to improve the accuracy of the code in space we “only” need to replace the function with the spatial differencing, and to improve the accuracy in time we “only” need to replace the solution of the ODE in time. This can also make it easier to test multiple methods rapidly.

Runge-Kutta

RK2:

\[ \begin{aligned} \symbf{\phi}^{[1]} &= \symbf{\phi}^{n} + \Delta t \, F ( \symbf{\phi}^n ) \, , \\ \symbf{\phi}^{n+1} &= \tfrac{1}{2} \left\{ \symbf{\phi}^{n} + \symbf{\phi}^{[1]} + \Delta t \, F ( \symbf{\phi}^{[1]} ) \right\} \, . \end{aligned} \]

from scipy.integrate import ode
def dphidt(t, phi):
    ...
r = ode(dphidt).set_integrator('dopri5', max_step=dt)
r.set_initial_value(...)
while t < t_end:
    r.integrate(r.t + dt)

The simplest Runge-Kutta update method is Euler, which is Forward Time differencing. The second order RK method here is one of a family of methods, but this choice has nice stability properties. Higher order RK methods exist, but get steadily less efficient above fourth order. This loss of efficiency tends to appear in two ways. First, more evaluations of the time-stepping function \(F\) are needed than you gain orders of accuracy (so four evaluations are needed for fourth order, but at least six are needed for fifth order). Second, the timestep restriction gets steadily worse with increasing order.

It is often most efficient to use library functions. Provided you can map your data to vectors at each stage (which should be straightforward with numpy operations or similar) then there should always be efficient, high accuracy library functions available, such as the scipy functions shown.

Summary

Semi-discretization simplifies constructing complex schemes.
Runge-Kutta methods are a standard approach
Libraries are an efficient starting point.