Jupyter notebooks

This is a Jupyter notebook using Python. You can install Jupyter locally to edit and interact with this notebook.

Finite difference methods for transient PDE

Method of Lines

Our method for solving time-dependent problems will be to discretize in space first, resulting in a system of ordinary differential equations

$$ M \dot u = f(u) $$

where the "mass matrix" $M$ might be diagonal and $f(u)$ represents a spatial discretization that has the form $f(u) = A u$ for linear problems.

Second Order Equations

A second order system such as

$$ \ddot u = f(t, u, \dot u) $$

may always be converted to a first order system by introducing a new variable

$$ \begin{bmatrix} u_0 \\ u_1 \end{bmatrix}' = \begin{bmatrix} u_1 \\ f(t, u_0, u_1) \end{bmatrix} . $$

Therefore, without loss of generality, we will focus on first order systems.

Implicit formulations

We have chosen the representation $M \dot u = f(t,u)$, but it is more general to write $h(t,u,\dot u) = 0$. If $\partial h/\partial \dot u$ is singular, then this describes a differential algebraic equation (DAE). General DAE are more challenging to solve and beyond the scope of this course, though we will allow a special case. We will consider cases in which the mass matrix $M$ is singular due to zero rows -- this is sometimes convenient to enforce boundary conditions and results in an "index 1" DAE.

Linear Equations

If $f(u) = A u$ is a linear function of $u$ then we have a linear ODE $$ \dot u = A(t) u + \text{source}(t) . $$ If $A(t)$ is independent of $t$ and $\text{source}(t) = 0$ then we have a linear, constant, autonomous ODE.

Non-autonomous ODE

Our ODE is "non-autonomous" if $\text{source}(t) \ne 0$. In that case, we can augment our system with a new independent variable $u_s$ and equation $u_s = 1$, then substitute $t = u_s$ to arrive at an equivalent autonomous formulation. The augmented equation will be integrated exactly using every method that we will consider. So there is no loss of generality by assuming that the system is autonomous, though many production ODE solvers will have direct support for non-autonomous systems.

Exponentials

When $u$ is a scalar then $A = \lambda$ is a scalar and the solution is $$ u(t) = u(0) e^{\lambda t} . $$

Question

What qualitative dynamics does this imply for

• $\lambda > 0$?
• $\lambda < 0$?
• $\lambda$ imaginary?

Question

• What if $A$ is diagonal?
• What if $A = X \Lambda X^{-1}$?

Matrix exponential

The general solution can be written in terms of the matrix exponential. $$ u(t) = e^{At} u(0) . $$ The matrix exponential is defined by its Taylor series $$ e^A = 1 + A + \frac{A^2}{2} + \frac{A^3}{3!} + \dotsb $$ and there are many practical ways to compute it.

Question

Suppose that the diagonalization $A = X \Lambda X^{-1}$ exists and derive a finite expression for the matrix exponential using the scalar exp function.

Forward Euler Method

The simplest method for solving $\dot u(t) = f(t,u)$ is $$ \tilde u(h) = u(0) + h f(0, u(0)) . $$

Let's try this on a scalar problem

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

where $k$ is a parameter controlling the rate at which the solution $u(t)$ is pulled toward the curve $\cos t$.

Questions

• What happens when $h$ is increased?
• What if $k$ is increased?

Example: Linear system

Now we consider linear systems

$$ \dot u = A u $$

which have an exact solution $u(t) = e^{At} u(0)$ in terms of the matrix exponential.

Questions

• Does shrinking $h$ make this more accurate?
• What if tfinal is extended?
• What are the eigenvalues of $A$?

Implicit methods

Recall that forward Euler is the step $$ \tilde u(h) = u(0) + h f(0, u(0)) . $$ This can be evaluated explicitly; all the terms on the right hand side are known so the approximation $\tilde u(h)$ is computed merely by evaluating the right hand side. Let's consider an alternative, backward Euler (or "implicit Euler"), $$ \tilde u(h) = u(0) + h f(h, \tilde u(h)) . $$ This is a (generally) nonlinear equation for $\tilde u(h)$.

We'll discuss the nonlinear case later and consider the linear case $f(t,u) = A u$, for which backward Euler reduces to

$$ \tilde u(h) = u(0) + h A \tilde u(h) $$

or

$$ (I - h A) u(h) = u(0) . $$

This requires solving a linear system involving both $h$ and $A$. Depending on the problem, this may be significantly more expensive than merely multiplying by $A$.

Midpoint Method

What if instead of evaluating the function at the end of the time step, we evaluated in the middle of the time step using the average of the endpoint values.

$$ \tilde u(h) = u(0) + h f\left(\frac h 2, \frac{\tilde u(h) + u(0)}{2} \right) $$

For the linear problem, this reduces to

$$ \Big(I - \frac h 2 A \Big) u(h) = \Big(I + \frac h 2 A\Big) u(0) .$$

Linear Stability Analysis

Why did forward Euler diverge while backward Euler decayed and Midpoint is apparently quite accurate for long time periods? We can answer this by considering the test equation $$ \dot u = \lambda u $$ and applying each method to construct $$ u(h) = R(h\lambda) u(0) $$ where $R(z)$ is called the stability function.

$$\begin{align} R(z) &= 1+z & \text{Forward Euler} \\ R(z) &= \frac{1}{1-z} & \text{Backward Euler} \\ R(z) &= \frac{1+z/2}{1-z/2} & \text{Midpoint} \end{align}$$

$\theta$ method

The above methods are all special cases of the $\theta$ method

$$ \tilde u(h) = u(0) + h f\left(\theta h, \theta\tilde u(h) + (1-\theta)u(0) \right) $$

which, for linear problems, is solved as

$$ (I - h \theta A) u(h) = \Big(I + h (1-\theta) A \Big) u(0) . $$

$\theta=0$ is explicit Euler, $\theta=1$ is implicit Euler, and $\theta=1/2$ is the midpoint rule. The stability function is $$ R(z) = \frac{1 + (1-\theta)z}{1 - \theta z}. $$

We will generalize slightly to allow solution of a linear differential algebraic equation

$$ M \dot u = A u + f(t,x) $$

where $M$ is (for now) a diagonal matrix that has zero rows at boundary conditions. With this generalization, the $\theta$ method becomes

$$ (M - h \theta A) u(h) = \Big(M + h (1-\theta) A \Big) u(0) + h f(h\theta, x) . $$

We will assume that $M$ is nonsingular if $\theta=0$.

Stiff decay to cosine

Observations

• $\theta=1$ is robust
• $\theta=1/2$ gets correct long-term behavior, but has oscillations at early times
• $\theta < 1/2$ allows oscillations to grow

Convergence plots

Observations

• $\theta=1$ is robust with large time steps (or large $k$).
• $\theta=1$ absolute error increases as initial transient is resolved.
• Only $\theta=1/2$ is second order accurate.
• $\theta < 1/2$ is unstable (diverges) for large time steps.

Definition: $A$-stability

A method is $A$-stable if the stability region $$ \{ z : |R(z)| \le 1 \} $$ contains the entire left half plane $$ \Re[z] \le 0 .$$ This means that the method can take arbitrarily large time steps without becoming unstable (diverging) for any problem that is indeed physically stable.

Definition: $L$-stability

A time integrator with stability function $R(z)$ is $L$-stable if $$ \lim_{z\to\infty} R(z) = 0 .$$ For the $\theta$ method, we have $$ \lim_{z\to \infty} \frac{1 + (1-\theta)z}{1 - \theta z} = \frac{1-\theta}{\theta} . $$ Evidently only $\theta=1$ is $L$-stable.

Oscillatory linear system

Observations

• If $h$ is too large, we don't resolve the solution.
• $\theta=1/2$ is the only second order method and is always more accurate than other choices, even with large time steps.

Transient PDE

Diffusion (heat equation)

Let's first consider diffusion of a quantity $u(t,x)$

$$ \dot u(t,x) - u''(t,x) = f(t,x) \qquad t > 0, -1 < x < 1 \\ u(0,x) = g(x) \qquad u(t,-1) = h_L(t) \qquad u'(t,1) = h_R(t) .$$

Let's use a Chebyshev discretization in space.

Observations

• Sharp central spike is diffused very quickly.
• Artifacts with $\theta < 1$.

Manufactured solution

Observations

• Errors are limited by time (not spatial) discretization error. This is a result of using the (spectrally accurate) Chebyshev method in space.
• $\theta=1$ is more accurate than $\theta = 1/2$, despite the latter being second order accurate in time. This is analogous to the stiff relaxation to cosine test.

Largest eigenvalues

Finite difference method

Question: max explicit Euler time step

Express the maximum stable time step $\Delta t$ using explicit Euler in terms of the grid spacing $\Delta x$.

Hyperbolic (wave) equations

The simplest hyperbolic equation is linear advection

$$ \dot u(t,x) + c u'(t,x) = f(t,x) $$

where $c$ is the wave speed and $f$ is a source term. In the homogenous ($f = 0$) case, the solution is given by characteristics

$$ u(t,x) = u(0, x - ct) . $$

This PDE also requires boundary conditions, but as a first-order equation, we can only enforce boundary conditions at one boundary. It turns out that this needs to be the inflow boundary, so if $c > 0$, that is the left boundary condition $u(t, -1) = g(t)$. We can solve this system using Chebyshev methods.

Observations

• $\theta > 1/2$ causes decay in amplitude
• $\theta < 1/2$ causes growth -- unstable
• An undershoot develops behind the traveling wave and increasing resolution doesn't make it go away
• We need an upwind boundary condition, otherwise the system is unstable
• Only Dirichlet inflow conditions are appropriate -- Neumann conditions produce a singular matrix

Finite difference

Observations

• Centered methods have an undershoot behind the traveling wave
• Upwind biasing of the stencil tends to reduce artifacts, but only stencil=2 removes undershoots
• Downwind biasing is usually unstable
• With upwinded stencil=2, we can use an explicit integrator, but the time step must satisfy $$ c \Delta t < \Delta x $$
• The upwind methods are in general dissipative -- amplitude is lost even with very accurate time integration
• The higher order upwind methods always produce artifacts for sharp transitions

Phase analysis

We can apply the advection differencing stencils to the test functions $$ \phi(x, \theta) = e^{i \theta x}$$ and compare to the exact derivative $$ \frac{d \phi}{d x} = i \theta \phi(x, \theta) . $$

Runge-Kutta methods

The methods we have considered thus far can all be expressed as Runge-Kutta methods, which are expressed in terms of $s$ "stage" equations (possibly coupled) and a completion formula. For the ODE

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

the Runge-Kutta method is

$$\begin{split} Y_i = u(t) + h \sum_j a_{ij} f(t+c_j h, Y_j) \\ u(t+h) = u(t) + h \sum_j b_j f(t+c_j h, Y_j) \end{split}$$

where $c$ is a vector of abscissa, $A$ is a table of coefficients, and $b$ is a vector of completion weights. These coefficients are typically expressed in a Butcher Table $$ \left[ \begin{array}{c|c} c & A \\ \hline & b^T \end{array} \right] = \left[ \begin{array}{c|cc} c_0 & a_{00} & a_{01} \\ c_1 & a_{10} & a_{11} \\ \hline & b_0 & b_1 \end{array} \right] . $$ We will see that, for consistency, the abscissa $c$ are always the row sums of $A$ and that $\sum_i b_i = 1$.

If the matrix $A$ is strictly lower triangular, then the method is explicit (does not require solving equations). We have seen forward Euler

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

backward Euler $$ \left[ \begin{array}{c|c} 1 & 1 \\ \hline & 1 \end{array} \right] ,$$ and Midpoint $$ \left[ \begin{array}{c|c} \frac 1 2 & \frac 1 2 \\ \hline & 1 \end{array} \right]. $$

Indeed, the $\theta$ method is $$ \left[ \begin{array}{c|c} \theta & \theta \\ \hline & 1 \end{array} \right] $$ and an alternative "endpoint" variant of $\theta$ (a generalization of the trapezoid rule) is $$ \left[ \begin{array}{c|cc} 0 & 0 & 0 \\ 1 & 1-\theta & \theta \\ \hline & 1-\theta & \theta \end{array} \right]. $$

Stability

To develop an algebraic expression for stability in terms of the Butcher Table, we consider the test equation

$$ \dot u = \lambda u $$

and apply the RK method to yield

$$ \begin{split} Y_i = u(0) + h \sum_j a_{ij} \lambda Y_j \\ u(h) = u(0) + h \sum_j b_j \lambda Y_j \end{split} $$

or, in matrix form,

$$ \begin{split} Y = \mathbb 1 u(0) + h \lambda A Y \\ u(h) = u(0) + h \lambda b^T Y \end{split} $$

where $\mathbb 1$ is a column vector of length $s$ consisting of all ones. This reduces to $$ u(h) = \underbrace{\Big( 1 + h\lambda b^T (I - h \lambda A)^{-1} \mathbb 1 \Big)}_{R(h\lambda)} u(0) . $$

Evidently the endpoint variant of $\theta$ has the same stability function as the original (midpoint) variant that we've been using. These methods are equivalent for linear problems, but different for nonlinear problems.

Higher order explicit methods: Heun's and RK4

Explicit Euler steps can be combined to create more accurate methods. One such example is Heun's method, $$ \left[ \begin{array}{c|cc} 0 & 0 & 0 \\ 1 & 1 & 0 \\ \hline & \frac 1 2 & \frac 1 2 \end{array} \right]. $$

Another explicit method is the famous four-stage RK4, $$ \left[ \begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0 \\ \frac 1 2 & \frac 1 2 & 0 & 0 & 0 \\ \frac 1 2 & 0 & \frac 1 2 & 0 & 0 \\ 1 & 0 & 0 & 1 & 0 \\ \hline & \frac 1 6 & \frac 1 3 & \frac 1 3 & \frac 1 6 \end{array} \right] . $$

Finally a method with lots of stability along the imaginary axis. Let's try it on some test problems.

Observations

• Solutions look pretty good and we didn't need a solve.
• We needed to evaluate the right hand side $s$ times per step

Work-precision diagrams for comparing methods

Since these methods do not cost the same per step, it is more enlightening to compare them using some measure of cost. For large systems of ODE, such as arise by discretizing a PDE, the cost of time integration is dominated by evaluating the right hand side (discrete spatial operator) on each stage. Measuring CPU time is a more holistic measure of cost, but the results depend on the implementation, computer, and possible operating system interference/variability. Counting right hand side function evaluations is a convenient, reproducible measure of cost.

Evidently Heun becomes resolved at lower cost than RK4.

Refinement in space and time

When solving a transient PDE, we should attempt to balance spatial discretization error with temporal discretization error. If we wish to use the same type of method across a range of accuracies, we need to

1. choose spatial and temporal discretizations with the same order of accuracy,

• choose grid/step sizes so the leading error terms are of comparable size, and
• ensure that both spatial and temporal discretizations are stable throughout the refinement range.

Since temporal discretization errors are proportional to the duration, simulations that run for a long time will need to use more accurate time discretizations.

Runge-Kutta order conditions

We consider the autonomous differential equation

$$ \dot u = f(u) . $$

Higher derivatives of the exact soultion can be computed using the chain rule, e.g.,

\begin{align*} \ddot u(t) &= f'(u) \dot u = f'(u) f(u) \\ \dddot u(t) &= f''(u) f(u) f(u) + f'(u) f'(u) f(u) . \\ \end{align*}

Note that if $f(u)$ is linear, $f''(u) = 0$. Meanwhile, the numerical solution is a function of the time step $h$,

$$\begin{split} Y_i(h) &= u(0) + h \sum_j a_{ij} f(Y_j) \\ U(h) &= u(0) + h \sum_j b_j f(Y_j). \end{split}$$

We will take the limit $h\to 0$ and equate derivatives of the numerical solution. First we differentiate the stage equations,

\begin{split} Y_i(0) &= u(0) \\ \dot Y_i(0) &= \sum_j a_{ij} f(Y_j) \\ \ddot Y_i(0) &= 2 \sum_j a_{ij} \dot f(Y_j) \\ &= 2 \sum_j a_{ij} f'(Y_j) \dot Y_j \\ &= 2\sum_{j,k} a_{ij} a_{jk} f'(Y_j) f(Y_k) \\ \dddot Y_i(0) &= 3 \sum_j a_{ij} \ddot f (Y_j) \\ &= 3 \sum_j a_{ij} \Big( \sum_k f''(Y_j) \dot Y_j \dot Y_k + f'(Y_j) \ddot Y_j \Big) \\ &= 3 \sum_{j,k,\ell} a_{ij} a_{jk} \Big( a_{j\ell} f''(Y_j) f(Y_k) f(Y_\ell) + 2 a_{k\ell} f'(Y_j) f'(Y_k) f(Y_\ell) \Big) \end{split}

where we have used Liebnitz's formula for the $m$th derivative, $$ (h \phi(h))^{(m)}|_{h=0} = m \phi^{(m-1)}(0) .$$ Similar formulas apply for $\dot U(0)$, $\ddot U(0)$, and $\dddot U(0)$, with $b_j$ in place of $a_{ij}$.

Equating terms $\dot u(0) = \dot U(0)$ yields $$ \sum_j b_j = 1, $$ equating $\ddot u(0) = \ddot U(0)$ yields $$ 2 \sum_{j,k} b_j a_{jk} = 1 , $$ and equating $\dddot u(0) = \dddot U(0)$ yields the two equations \begin{split} 3\sum_{j,k,\ell} b_j a_{jk} a_{j\ell} &= 1 \\ 6 \sum_{j,k,\ell} b_j a_{jk} a_{k\ell} &= 1 . \end{split}

Observations

• These are systems of nonlinear equations for the coefficients $a_{ij}$ and $b_j$. There is no guarantee that they have solutions.
• The number of equations grows rapidly as the order increases.

	$u^{(1)}$	$u^{(2)}$	$u^{(3)}$	$u^{(4)}$	$u^{(5)}$	$u^{(6)}$	$u^{(7)}$	$u^{(8)}$	$u^{(9)}$	$u^{(10)}$
# terms	1	1	2	4	9	20	48	115	286	719
cumulative	1	2	4	8	17	37	85	200	486	1205

• Usually the number of order conditions does not exactly match the number of free parameters, meaning that the remaining parameters can be optimized (usually numerically) for different purposes, such as to minimize the leading error terms or to maximize stability in certain regions of the complex plane. Finding globally optimal solutions can be extremely demanding.
• The arithmetic managing the derivatives gets messy, but can be managed using rooted trees.

Theorem (from Hairer, Nørsett, and Wanner)

A Runge-Kutta method is of order $p$ if and only if $$ \gamma(\mathcal t) \sum_{j} b_j \Phi_j(t) = 1 $$ for all trees $t$ of order $\le p$.

For a linear autonomous equation $$ \dot u = A u $$ we only need one additional order condition per order of accuracy because $f'' = 0$. These conditions can also be derived by equating derivatives of the stability function $R(z)$ with the exponential $e^z$. For a linear non-autonomous equation $$ \dot u = A(t) u + g(t) $$ or more generally, an autonomous system with quadratic right hand side, $$ \dot u = B (u \otimes u) + A u + C $$ where $B$ is a rank 3 tensor, we have $f''' = 0$, thus limiting the number of order conditions.

Embedded error estimation and adaptive control

It is often possible to design Runge-Kutta methods with multiple completion orders, say of order $p$ and $p-1$.

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

The classical RK4 does not come with an embedded method, but most subsequent RK methods do.

The Bogacki-Shampine method is given by

While this method has four stages, it has the "first same as last" (FSAL) property meaning that the fourth stage exactly matches the completion formula, thus the first stage of the next time step. This means it can be implemented using only three function evaluations per time step.

Higher order methods with embedded error estimation include

• Fehlberg, a 6-stage, 5th order method for which the 4th order embedded formula has been optimized for accuracy.
• Dormand-Prince, a 7-stage, 5th order method with the FSAL property, with the 5th order completion formula optimized for accuracy.

Adaptive control

Given a completion formula $b^T$ of order $p$ and $\tilde b^T$ of order $p-1$, an estimate of the local truncation error (on this step) is given by $$ e_{\text{loc}}(h) = \lVert h (b - \tilde b)^T f(Y) \rVert \in O(h^p) . $$ Given a tolerance $\epsilon$, we would like to find $h_*$ such that $$ e_{\text{loc}}(h_*) < \epsilon . $$ If $$e_{\text{loc}}(h) = c h^p$$ for some constant $c$, then $$ c h_*^p < \epsilon $$ implies $$ h_* < \left( \frac{\epsilon}{c} \right)^{1/p} . $$ Given the estimate with the current $h$, $$ c = e_{\text{loc}}(h) / h^p $$ we conclude $$ \frac{h_*}{h} < \left( \frac{\epsilon}{e_{\text{loc}}(h)} \right)^{1/p} . $$

Notes

• Usually a "safety factor" less than 1 is included so the predicted error is less than the threshold to reject a time step.
• We have used an absolute tolerance above. If the values of solution variables vary greatly in time, a relative tolerance $e_{\text{loc}}(h) / \lVert u(t) \rVert$ or a combination thereof is desirable.
• There is a debate about whether one should optimize the rate at which error is accumulated with respect to work (estimate above) or with respect to simulated time (as above, but with error behaving as $O(h^{p-1})$). For problems with a range of time scales at different periods, this is usually done with respect to work.
• Global error control is an active research area.

Homework 2: Due 2017-10-18 (Wednesday)

• Implement an explicit Runge-Kutta integrator that takes an initial time step $h_0$ and an error tolerance $\epsilon$.
• You can use the Bogacki-Shampine method or any other method with an embedded error estimate.
• A step should be rejected if the local truncation error exceeds the tolerance.
• Test your method on the nonlinear equation $$ \begin{bmatrix} \dot u_0 \\ \dot u_1 \end{bmatrix} = \begin{bmatrix} u_1 \\ k (1-u_0^2) u_1 - u_0 \end{bmatrix} $$ for $k=2$, $k=5$, and $k=20$.
• Make a work-precision diagram for your adaptive method and for constant step sizes.
• State your conclusions or ideas (in a README, or Jupyter notebook) about appropriate (efficient, accurate, reliable) methods for this type of problem.

Implicit Runge-Kutta methods

We have been considering examples of high-order explicit Runge-Kutta methods. For processes like diffusion, the time step becomes limited (under grid refinement, but usually for practical resolution) by stability rather than accuracy. Implicit methods, especially $A$-stable and $L$-stable methods, allow much larger time steps.

Diagonally implicit

A Runge-Kutta method is called diagonally implicit if the Butcher matrix $A$ is lower triangular, in which case the stages can be solved sequentially. Each stage equation has the form $$ Y_i - h a_{ii} f(Y_i) = u(0) + h \sum_{j<i} a_{ij} f(Y_j) $$ where all terms in the right hand side are known. For stiff problems, it is common to multiply though by $\alpha = (h a_{ii})^{-1}$, yielding $$ \alpha Y_i - f(Y_i) = \alpha u(0) + \sum_{j<i} \frac{a_{ij}}{a_{ii}} f(Y_j) . $$

• It is common for solvers to reuse a linearization associated with $f(Y_i)$.
• It is common to have setup costs associated with the solution of the "shifted" problem.

Methods with constant diagonals, $a_{ii} = a_{jj}$, are often desired to amortize setup costs. These methods are called singly diagonally implicit. There are also related methods called Rosenbrock or Roserbrock-W that more aggressively amortize setup costs.