44  A6: Solutions

using LinearAlgebra, Plots, LaTeXStrings

In lectures, we have looked at

\begin{align} \tag{IVP} &u'(t) = f(t, u(t)) \nonumber\\ \tag{BVP} &-u'' + p u' + qu = f\nonumber\\ \tag{heat} &u_t - \alpha^2 \Delta u = f \end{align}

and various methods to approximate the solution (IVP: Euler, one-step, Runge-Kutta, multistep; BVP: finite differences (including with FFT), finite elements; Heat equation: finite differences).

In this assignment, you’ll analyse simple finite difference approximations to the linear transport equation:

\begin{align} &u_t + c u_x = 0 \nonumber\\ &u(x,0) = u_0(x) \tag{transport} \end{align}

with constant velocity c \in \mathbb R. In lectures, we saw that the solution is given by u(x,t) = u_0(x-ct). While in this case the exact solution is known, it is still useful to analyse methods for approximating this solution numerically. These methods apply more generally to more complicated hyperbolic equations.

Consider the spatial domain [a,b] and temporal domain [0,T] with meshes given as follows: x_i := a + \frac{b-a}{m} i for i = 0,\dots,m and t_j := j\frac{T}{n} for j = 0,\dots,n. Again, we define h := \frac{b-a}{m} and \tau := \frac{T}{n}. We aim to approximate u(x_i, t_j) with u_{i,j}.

Exercise 1. Suppose one approximates the linear transport equation with simple forward differences for both the space and time derivatives. Show that this approximation gives

\begin{align} u_{i,j+1} = (1+\sigma) u_{i,j} - \sigma u_{i+1,j} \tag{transport$_{\tau, h}$} \end{align}

where \sigma:= \frac{c \tau}{h}.

We approximate

\begin{align} u_t(x_i,t_j) &\approx \frac{u_{i,j+1} - u_{i,j}}{\tau} \nonumber\\ u_{x}(x_i,t_j) &\approx \frac{u_{i+1,j} - u_{i,j}}{h} \end{align}

and thus we approximate the PDE u_t + c u_x = 0 with

\begin{align} \frac{u_{i,j+1} - u_{i,j}}{\tau} + c \frac{u_{i+1,j} - u_{i,j}}{h} &= 0. \end{align}

Multiplying by \tau gives

\begin{align} u_{i,j+1} - u_{i,j}+ \frac{c\tau}{h}( u_{i+1,j} - u_{i,j} ) &= 0, \end{align}

and thus we have \text{(transport}_{\tau,h}\text{)}.

Remark (Consistency).

Notice that, if u is a sufficiently smooth solution to u_t + c u_x = 0, then

\begin{align} &\frac{u(x_i,t_{j+1}) - u(x_i,t_j)}{\tau} + c \frac{u(x_{i+1},t_j) - u(x_i,t_j)}{h} \nonumber\\ % &= \left( u_t + \frac{\tau}{2} u_{tt} + \frac{\tau^2}{6} u_{ttt} + ... \right) + c \left( u_x + \frac{h}{2} u_{xx} + \frac{h^2}{6} u_{xxx} + ... \right) \nonumber\\ % &= \left( u_t + c u_x \right) + \frac{1}{2} \left( \tau u_{tt} + c h u_{xx} \right) + \frac{1}{6} \left( \tau^2 u_{ttt} + ch^2 u_{xxx} \right) % + \cdots \end{align}

Here, the functions are evaluated at (x_i,t_{j}). Therefore, as h,\tau\to0, the local truncation error goes to zero at least as quickly as O(\tau) + O(h).

When u is a solution to u_t + c u_x = 0, then we have u_{tt} = (u_t)_t % = (- c u_x)_t % = -c (u_t)_x % = -c ( -c u_x )_x % = c^2 u_{xx} and so the leading order term is

\begin{align} \frac{1}{2} ( \tau u_{tt} + c h u_{xx} ) &= \frac{c h}{2} \left(\sigma + 1\right) u_{xx} \end{align}

(where \sigma = \frac{c\tau}{h}). As a result we expect that, if the method is stable, then the order of convergence is O(h) in general and O(h^2) if \sigma = -1.

Exercise 2. We consider the linear transport equation with periodic boundary conditions: u(a,t) = u(b,t) and u_x(a,t) = u_x(b,t). Complete the following code sample to implement the finite difference scheme from exercise 1:

function fdm_transport(m, n, c, u0; a=0., b=1., T=5, warn=true)
    h, τ = (b-a)/m, T/n 
    x, t = collect(a:h:b), collect(0:τ:T)
    
    u = zeros(m+1, n+1)
    # Initial condition
    u[:, 1] = u0.(x)

    # CLF
    σ = c * τ / h 
    if warn 
        println("σ = $(round(σ, digits=2))")
    end

    for j ∈ 1:n 
        for i ∈ 1:m 
            # periodic neighbour of m is 1:
            i′ = (i%m) + 1

            # finite difference rule here:
            u[i,j+1] = (1+σ)u[i,j] - σ*u[i′,j] 
        end

        # periodic boundary
        u[m+1,j+1] = u[1,j+1]
    end

    return x, t, u
end

a, b, T = -10, 10, 20
m, n, c = 100, 100, -1.

# initial condition (w/ PBCs)
g(x) = sin(x) * exp(-x^2)
u0(x) = g(a + mod(x - a, b-a))

# exact solution 
u_e(x, t) = u0(x - c*t)

x,t,u = fdm_transport(m, n, c, u0; a, b, T)

anim = @animate for j ∈ 1:length(t)
    scatter(x, u[:, j]; 
        title = L"Transport equation $u_t + c u_{x} = 0$",
        label = "t = $(t[j])", legend=:topright,
        ylims = [-1/2,1/2], m=2)
    plot!( x->u_e(x,t[j]); 
        label="exact solution", linestyle=:dot)
end
gif(anim, "pics/transport-1.gif", fps=10)

σ = -1.0

[ Info: Saved animation to c:\Users\math5\Math 5485\Math5486\pics\transport-1.gif

Exercise 3. Experiment with the above code. What happens when

• (i) \sigma < -1,
  
• (ii) \sigma = -1,
  
• (iii) -1 < \sigma < 0,
  
• (iv) \sigma = 0, and
  
• (v) \sigma>0.

Either plot the numerical solution in each case, or write down the parameters you used.

For all choices of parameters that I tested, we see the following behaviour:

• \sigma < -1 : unstable,
  
• \sigma = -1 : stable and the numerical solution agrees with the exact solution for long times (see graph above),
  
• \sigma \in (-1,0) : stable but the solution decays as a function of time (this is known as numerical diffusion or numerical viscosity; the numerical scheme introduces a “ghost” diffusive term - see below remark). When \sigma is closer to zero, this decay becomes more noticable,
  
• \sigma = 0 : stable, the solution and numerical solution are both stationary,
  
• \sigma > 0 : unstable.

σ = -1.1

[ Info: Saved animation to c:\Users\math5\Math 5485\Math5486\pics\transport-2.gif

σ = -0.9

[ Info: Saved animation to c:\Users\math5\Math 5485\Math5486\pics\transport-3.gif

σ = -0.2

[ Info: Saved animation to c:\Users\math5\Math 5485\Math5486\pics\transport-4.gif

σ = 0.0

[ Info: Saved animation to c:\Users\math5\Math 5485\Math5486\pics\transport-5.gif

σ = 0.2

[ Info: Saved animation to c:\Users\math5\Math 5485\Math5486\pics\transport-6.gif

σ = 1.0

[ Info: Saved animation to c:\Users\math5\Math 5485\Math5486\pics\transport-7.gif

Remark (Numerical diffusion).

Recall from the “consistency” remark that, if \sigma = \frac{c \tau}{h} is fixed, then

\begin{align} &\frac{u(x_i,t_{j+1}) - u(x_i,t_j)}{\tau} + c \frac{u(x_{i+1},t_j) - u(x_i,t_j)}{h} \nonumber\\ % &= (u_t + c u_x) + \frac{c h}{2} (\sigma + 1) u_{xx} % + O(h^2) \end{align}

That is, our finite difference scheme is actually approximating the solution of the equation u_t + c u_x + \nu u_{xx} = 0 where \nu := \frac{c h}{2} (\sigma + 1). We have seen (in the heat equation) that \nu u_{xx} with \nu < 0 is a “diffusive” term that smooths out solutions and causes solutions to decay in time. Similarly, when \nu >0, this introduces instability (peaks in the solution have u_{xx}<0 which makes u increase in time!). The \nu u_{xx} term is known as numerical diffusion because it doesn’t come from the exact equation, it only comes from the numerical scheme used to approximate it.

Here, we explain the findings above:

• \sigma < -1 : unstable; we have c < 0 and so \nu > 0,
• \sigma = -1 : stable; no numerical diffusion,
  
• \sigma \in (-1,0) : \nu < 0 so the solution decays. When \sigma is closer to 0, \nu is more negative so the solution decays faster,
• \sigma = 0 : stable, the solution and numerical solution are both stationary,
  
• \sigma > 0 : unstable; \nu > 0.

Here are some error curves in the max-norm at time T with fixed \sigma = \frac{c\tau}{h}:

The domain of dependence of the solution at (x,t) is the set of all y for which changing u_0(y) can affect the solution u at (x,t).

Exercise 4. Write down the domain of dependence of the solution to the linear transport equation at (x,t).

The domain of dependence of the solution at (x,t) is \{ x - c t \} (since the solution is u(x,t) = u_0(x-ct)).

The numerical domain of dependence at (x_i,t_j) is the set of all \{x_k\} such that the initial data u_{k,0} can affect u_{i,j}.

Exercise 5. Write down the numerical domain of dependence of \text{(transport}_{h,\tau}\text{)}.

Since u_{i,j} = (1+\sigma) u_{i,j-1} - \sigma u_{i+1,j-1}, we find that u_{ij} depends on u_{i,j-1} and u_{i+1,j-1}. Therefore, u_{i,j-1} depends on u_{i,j-2} and u_{i+1,j-2} and u_{i+1,j-1} depends on u_{i+1,j-2} and u_{i+2,j-2}. Iterating this process, we find that u_{ij} depends on u_{i,j-p},u_{i+1,j-p}, \cdots, u_{i+p,j-p} for all p\in\mathbb N. As a result, the numerical domain of dependence is \{x_{i},x_{i+1},\dots,x_{i+j}\}.

The Courant–Friedrichs–Lewy (CFL) condition states that, for a numerical method to converge, the exact domain of dependence must be contained in the numerical domain of dependence as h,\tau \to 0.

Exercise 6. Under what condition on \sigma does our numerical method for the linear transport equation satisfy the CFL condition.

The CFL condition reads x_i \leq x_i - c t_j \leq x_{i+j} which is equivalent to \sigma \in [-1,0] where \sigma := \frac{c \tau}{h}.

Notice. This is the same condition as we saw before!

Remark: If the CLF condition is violated, then one can change some of the initial data that is both outside the domain of numerical dependence (so does not affect the numerical solution) and inside the domain of exact dependence (so does affect the exact solution). As a result, the numerical scheme cannot approximate the true solution. The CLF condition is necessary but not sufficient for a numerical scheme to converge (as we will now see):

Consider the following centered difference approximation to the linear transport equation:

\begin{align} u_{i,j+1} = -\frac12 \sigma u_{i+1,j} + u_{ij} + \frac{1}2\sigma u_{i-1,j} \end{align}

where \sigma = \frac{c\tau}{h}.

Exercise 7. Show that, for |\sigma| \leq 1, this scheme satisfies the CFL condition.

In this case the domain of numerical dependence is \{x_{i-j},\dots,x_{i+j}\} and so the CFL condition is

\begin{align} x_{i-j} \leq x_i - c t_j \leq x_{i+j} \end{align}

which is equivalent to x_i - j h \leq x_i - c j \tau \leq x_i + j h and $ - 1 $.

Exercise 8. Show that this scheme is unstable for all \sigma\not=0 via a von Neumann stability analysis.

We take u_{ij} := e^{\mathrm{i} k x_i} and substitute this into the equation:

\begin{align} u_{i,j+1} &= -\frac12 \sigma e^{\mathrm{i} k x_{i+1}} + e^{\mathrm{i} k x_i} + \frac{1}2\sigma e^{\mathrm{i} k x_{i-1}} \nonumber\\ % &= \left( 1 -\frac12 \sigma e^{\mathrm{i} k h} + \frac{1}2\sigma e^{-\mathrm{i} k h } \right) e^{\mathrm{i} k x_{i}} \nonumber\\ % &= \left( 1 -\frac12 \sigma \left[ e^{\mathrm{i} k h} - e^{-\mathrm{i} k h } \right] \right) e^{\mathrm{i} k x_{i}} \nonumber\\ % &= \left( 1 - \mathrm{i} \sigma \sin(kh) \right) e^{\mathrm{i} k x_{i}} \nonumber\\ % &=: \lambda(k) e^{\mathrm{i} k x_i} . \end{align}

Here, we have

\begin{align} |\lambda(k)|^2 = 1 + \sigma^2 \sin^2(kh) > 1 \end{align}

for all k \not\in \frac{2\pi}{h} \mathbb Z. As a result, this method is unconditionally unstable.

Here, we implement this method (and we can see that it is unstable for our choice of parameters):

function fdm_transport_unstable(m, n, c, u0; a=0., b=1., T=5)
    h, τ = (b-a)/m, T/n 
    x, t = collect(a:h:b), collect(0:τ:T)
    
    u = zeros(m+1, n+1)
    u[:, 1] = u0.(x)

    # CLF
    σ = c * τ / h 
    println("σ = $(round(σ, digits=2))")

    for j ∈ 1:n 
        for i ∈ 1:m 
            i₋ = (i==1) ? m : i-1
            i₊ = (i==m) ? 1 : i+1
            u[i,j+1] = (1/2)*σ*u[i₋,j] + u[i,j] + (-1/2)*σ*u[i₊,j] 
        end
        u[m+1,j+1] = u[1,j+1]
    end

    return x, t, u
end

a, b, T = -10, 10, 20
m, n, c = 100, 100, .5

# initial condition (w/ PBCs)
g(x) = sin(x) * exp(-x^2)
u0(x) = g(a + mod(x - a, b-a))

# exact solution 
u_e(x, t) = u0(x - c*t)

x,t,u = fdm_transport_unstable(m, n, c, u0; a, b, T)

anim = @animate for j ∈ 1:length(t)
    plot(x, u[:, j]; 
        title = L"Transport equation $u_t + c u_{x} = 0$ (unstable)",
        label = "t = $(t[j])", legend=:topright,
        ylims = [-1/2,1/2], m=2)
    plot!( x->u_e(x,t[j]); 
        label="exact solution", linestyle=:dot)
end
gif(anim, "pics/transport-2.gif", fps=10)

σ = 0.5

[ Info: Saved animation to c:\Users\math5\Math 5485\Math5486\pics\transport-2.gif