# Simple time-stepping#

Our previous examples were all stationary problems. However, many practical simulations describe processes that change over time. In this part we want to start looking a bit closer onto time-domain partial differential equations and their efficient implementation.

## Finite Difference Approximation for the time-derivative#

We want to approximate the derivative \(\frac{du}{dt}\). Remember that

for sufficiently small \(h\).

There are three standard approximations for the time-derivative:

The forward difference:

The backward difference:

The centered difference:

To understand the error of these schemes we can use Tayler expansions to obtain

Hence, the error of the first two schemes decreases linearly with \(h\) and the error in the centred scheme decreases quadratically with \(h\).

## The 3-point stencil for the second derivative#

For simplicity we denote \(u_i := u(t)\), \(u_{i + 1} := u(t + h)\), \(u_{i-1} := u(t - h)\). We want to approximate

The trick is to use an approximation around half-steps for the outer derivative, resulting in

The derivatives at the half-steps are now again approximated by centered differences, resulting in

This is the famous second order finite difference operator that we have already used before. Its error is quadratically small in \(h\).

## Application to time-dependent problems#

We now want to solve

where \(U(t):\mathbb{R}\rightarrow \mathbb{R}^n\) is some vector valued function.

The idea is replace \(\frac{dU}{dt}\) by a finite difference approximation.

Forward Euler Method

Backward Euler Method

The forward Euler method is an explicit method. We have that

and the right-hand side only has known values.

In contrast to this is the backward Euler Method, which is an implicit method since

We hence need to solve a linear or nonlinear system of equations to compute \(U_{n+1}\).

## Stability of forward Euler#

We consider the model problem

for \(\alpha < 0\). Note that the explicit solution of this problem is \(u(t) = u_0e^{\alpha t}\). For \(t\rightarrow\infty\) we have \(u(t)\rightarrow 0\) if \(\alpha < 0\).

The forward Euler method can now be written as

Hence, in order for the solution to decay we need that \(|1+\alpha\Delta t| < 1\) or equivalently

from which we obtain \(|\alpha\Delta t| < 2\) (if \(\alpha\) negative). Now consider the problem

with \(A\in\mathbb{R}^{n\times n}\) diagonalizable. For any eigenpair \((\lambda, \hat{U})\) of \(A\) satisfying \(A\hat{U} = \lambda\hat{U}\) the function \(U(t) = e^{\lambda t}\hat{U}\) is a solution for this problem. If the eigenvalues are real and negative we require for forward Euler to converge that

where \(\lambda_{max}\) is the largest eigenvalue by magnitude. Note that if the eigenvalues are complex the condition becomes

where \(\lambda\) is the eigenvalue with largest negative real part by magnitude.

As example let us take a look at the problem

with \(u(x, 0) = u_0(x)\), \(u(0, t) = u(1, t) = 0\). We can discretise the right-hand side using our usual second order finite dfference scheme. For the left-hand side, we use the forward Euler method. This gives us the recurrence equation

with \(A = \frac{1}{h^2}\text{tridiag}(1, -2, 1)\).

The eigenvalues of \(A\) are given explicitly as

We therefore have that \(|\lambda_{max}|\sim \frac{4}{h^2}\). Hence, for forward Euler to be stable we require that

Hence, we need that \(\Delta t\sim h^2\), meaning that the number of required time steps grows qudratically with the discretisation accuracy.

## Stability of backward Euler#

For backward Euler we obtain

We now require that \(|(1-\alpha \Delta t)^{-1}| < 1\). But for \(\alpha<0\) this is always true. Hence, the backward Euler method is unconditionally stable.

## Implicit vs explicit methods#

This analysis is very typical. In computational sciences we always have to make a choice between implicit and explicit methods. The advantage of implicit methods are the very good stability properties, allowing us for the backward Euler method to choose the time-discretisation independent of the spatial discretisation. For explicit methods we have to be much more careful and in the case of Euler we have the quadratic dependency between time-steps and spatial discretisation. However, a single time-step is much cheaper for explicit Euler as we do not need to solve a linear or nonlinear system of equations in each step. The right choice of solver depends on a huge number of factors and is very application dependent.

## Time-Stepping Methods in Software#

In practice we do not necesarily use explicit or implicit Euler. There are many better methods out there. The Scipy library provides a number of time-stepping algorithms. For PDE problems PETSc has an excellent infrastructure of time-stepping methods built in to support the solution of time-dependent PDEs.