Numerical PDE I: Finite Difference

Introduction¶

Partial Differential Equations (PDEs) are fundamental tools in the mathematical modeling of various physical phenomena. Unlike Ordinary Differential Equations (ODEs), which involve functions of a single variable and their derivatives, PDEs involve functions of multiple variables and their partial derivatives. This distinction makes PDEs particularly powerful in describing systems where changes occur in more than one dimension, such as in space and time.

A PDE is an equation that relates the partial derivatives of a multivariable function. In general form, a PDE can be written as:

\begin{align} F\left(x_1, x_2, \ldots, x_n, u, \frac{\partial u}{\partial x_1}, \frac{\partial u}{\partial x_2}, \ldots, \frac{\partial^k u}{\partial x_1^{k_1}\partial x_2^{k_2} \ldots \partial x_n^{k_n}}\right) = 0 \end{align}

(1)

where $u = u(x_1, x_2, \ldots, x_n)$ is the unknown function, and $\partial u/\partial x_i$ denotes the partial derivatives of $u$ with respect to the variables $x_i$ .

PDEs are essential in modeling continuous systems where the state of the system depends on multiple variables. They appear in various fields such as physics, engineering, finance, and biology, describing phenomena like heat conduction, wave propagation, fluid dynamics, and quantum mechanics.

Definition and Significance in Modeling Continuous Systems¶

PDEs provide a framework for formulating problems involving functions of several variables and their rates of change. They are indispensable in describing the behavior of physical systems where spatial and temporal variations are intrinsic. Examples include:

Advection Equation: Models a quantity, e.g., density, moves with velocity $c$ :
$\begin{align} \frac{\partial u}{\partial t} + c \nabla u = 0. \end{align}$
(2)
Heat Equation: Models the distribution of heat (or temperature) in a given region over time.
$\begin{align} \frac{\partial u}{\partial t} = \alpha \nabla^2 u. \end{align}$
(3)
Wave Equation: Describes the propagation of waves, such as sound or electromagnetic waves, through a medium.
$\begin{align} \frac{\partial^2 u}{\partial t^2} = c^2 \nabla^2 u \end{align}$
(4)
Laplace’s Equation: Represents steady-state solutions where the system does not change over time, such as electric potential in a region devoid of charge.
$\begin{align} \nabla^2 u = 0 \end{align}$
(5)

The ability to model such diverse phenomena underscores the versatility and importance of PDEs in scientific and engineering disciplines.

Finite Difference Methods (FDM)¶

Solving PDEs numerically is essential for modeling complex physical systems that lack closed-form analytical solutions. Among the various numerical methods available, Finite Difference Methods (FDM) are particularly popular due to their simplicity and ease of implementation. They approximate the derivatives in PDEs by using difference on a discretized grid. This approach transforms continuous PDEs into discrete algebraic equations that can be solved iteratively. FDM is widely used in engineering and scientific computations due to its straightforward application to regular grids and its compatibility with existing numerical solvers.

Forward Time Centered Space¶

The Forward Time Centered Space (FTCS) scheme is one of the simplest explicit finite difference methods used to solve time-dependent PDEs. It approximates the time derivative using a forward difference and the spatial derivatives using centered differences.

Consider the linear advection equation that models the transport of a quantity $u$ with constant speed $c$ :

\begin{align} \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0 \end{align}

(6)

Let $u_i^n$ denote the numerical approximation of $u$ at spatial index $i$ and time level $n$ . The FTCS scheme discretizes the advection equation as follows:

\begin{align} \frac{u_i^{n+1} - u_i^n}{\Delta t} + c \frac{u_{i+1}^n - u_{i-1}^n}{2 \Delta x} = 0 \end{align}

(7)

Solving for $u_i^{n+1}$ :

\begin{align} u_i^{n+1} = u_i^n - \frac{c \Delta t}{2 \Delta x} \left( u_{i+1}^n - u_{i-1}^n \right) \end{align}

(8)

This explicit update rule allows the computation of the solution at the next time step based on the current and neighboring spatial points.

Here is a simple python implementation:

# Parameters

c  = 1.0    # advection speed
l  = 1.0    # domain size
dt = 0.005  # time step

nx = 100    # number of spatial points
nt = 200    # number of time steps

import numpy as np

X, dx = np.linspace(0, l, nx, endpoint=False, retstep=True)  # spatial grid
U0    = np.sin(2*np.pi * X)  # initial condition

# Forward Time Centered Space (FTCS) scheme

def FTCS(c, U0, dx, dt, n):
    sigma = c * dt / dx
    U = [U0]
    for _ in range(n):
        U0 = U[-1]
        U1 = U0 - sigma * (np.roll(U0,-1) - np.roll(U0,1)) / 2
        U.append(U1)
    return np.array(U)

U     = np.sin(2*np.pi * (X - c*dt*nt)) # analytical solution
UFTCS = FTCS(c, U0, dx, dt, nt)         # numerical solution

Let’s now plot the result! After $t = dt n_t = 1$ , the solution should match the initial condition exactly.

from matplotlib import pyplot as plt

plt.plot(X, U0,              label='Initial Condition')
plt.plot(X, U, ':',          label='Exact Solution')
plt.plot(X, UFTCS[-1], '.-', label='FTCS Scheme')
plt.xlabel('x')
plt.ylabel('u')
plt.legend()

However, the numerical solution is oscillating. This looks like a numerical artifact. Let’s inspect it with a movie.

from matplotlib.animation import ArtistAnimation
from IPython.display import HTML
from tqdm import tqdm

def animate(X, U):
    fig, ax = plt.subplots(1,1)
    ax.set_xlabel('x')
    ax.set_ylabel('y')

    frames = []
    for n in tqdm(range(len(U))):
        f = ax.plot(X, U[n], 'C0.-', animated=True)
        frames.append(f)
        plt.close()
    
    return ArtistAnimation(fig, frames, interval=50)

anim = animate(X, UFTCS)

HTML(anim.to_html5_video())  # display animation
# anim.save('FTCS.mp4')        # save animation

The oscillation grows slowly as the solution evolve. This looks like numerical instability that we discussed in ODE integrator.

# HANDSON: change `c`, `dt`, `dx`, or other parameters to observe how
#          the numerical solution changes.

# HANDSON: change the initial condition to, e.g., a boxcar function,
#          and observe how the numerical solution changes

To understand the oscillation, we will perform a stability analysis.

Von Neumann Stability Analysis¶

Von Neumann stability analysis is a mathematical technique used to assess the stability of finite difference schemes applied to (linear) PDEs. This method involves decomposing the numerical solution into Fourier modes and analyzing the growth or decay of these modes over time. If all Fourier modes remain bounded, the numerical scheme is considered stable. Otherwise, the scheme is unstable and may produce erroneous results.

To perform von Neumann stability analysis, we assume a solution of the form:

\begin{align} u_i^n = G^n e^{ikx_i} \end{align}

(9)

where:

$u_i^n$ is the numerical approximation of $u$ at spatial index $i$ and time level $n$ .
$G$ is the amplification factor (recall our ODE stability analysis).
$k$ is the wave number.
$x_i = i \Delta x$ is the spatial position of the $i$ -th grid point.

The goal is to determine whether the magnitude of the amplification factor $|G|$ remains less than or equal to 1 for all possible wave numbers $k$ . If $|G| \leq 1$ for all $k$ , the numerical scheme is stable.

Now, consider the FTCS scheme applied to the linear advection equation

\begin{align} \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0, \end{align}

(10)

where the FTCS update rule reads

\begin{align} u_i^{n+1} = u_i^n - \frac{c\Delta t}{2\Delta x} \left(u_{i+1}^n - u_{i-1}^n\right). \end{align}

(11)

Assumed a single Fourier mode solution,

\begin{align} u_i^n = G^n e^{ikx_i}. \end{align}

(12)

Substituting it into the FTCS update rule,

\begin{align} G e^{ikx_i} = e^{ikx_i} - \frac{c\Delta t}{2\Delta x} \left(e^{ikx_{i+1}} - e^{ikx_{i-1}}\right). \end{align}

(13)

Simplify using $x_{i \pm 1} = x_i \pm \Delta x$ , we have

\begin{align} G = 1 - \frac{c \Delta t}{2\Delta x} \left(e^{ik\Delta x} - e^{-ik\Delta x}\right). \end{align}

(14)

Using Euler’s formula $e^{i\theta} - e^{-i\theta} = 2i \sin\theta$ :

\begin{align} G = 1 - i \frac{c \Delta t}{\Delta x} \sin(k \Delta x) = 1 - i \sigma \sin(k \Delta x) \end{align}

(15)

where $\sigma \equiv c \Delta t/\Delta x$ is the Courant number.

Calculating the magnitude of $G$ , we obtain:

\begin{align} |G|^2 = 1 + \sigma^2 \sin^2(k \Delta x) > 1 \end{align}

(16)

Since $|G| > 1$ for any $\sigma$ , the FTCS scheme is unconditionally unstable for the linear advection equation.

To overcome this limitation, we will introduce two more robust finite difference methods: the Upwind Scheme and the Lax-Wendroff Scheme. These methods enhance stability and accuracy, making them more suitable for solving advection-dominated problems.

Upwind Scheme¶

The Upwind Scheme is a finite difference method specifically designed to handle advection-dominated problems more effectively than symmetric schemes like FTCS. By incorporating the direction of wave propagation into the discretization of spatial derivatives, the upwind method enhances numerical stability and reduces non-physical oscillations.

In advection processes, information propagates in a specific direction determined by the flow velocity $c$ . The upwind scheme leverages this directional information to bias the spatial derivative approximation, ensuring that the numerical flux aligns with the physical transport direction. This directional bias significantly improves the stability of the numerical solution, especially when dealing with sharp gradients or discontinuities.

The upwind scheme discretizes the spatial derivative based on the sign of the advection speed $c$ :

For $c > 0$ (flow to the right):
$\begin{align} \frac{\partial u}{\partial x} \approx \frac{u_i^n - u_{i-1}^n}{\Delta x} \end{align}$
(17)
For $c < 0$ (flow to the left):
$\begin{align} \frac{\partial u}{\partial x} \approx \frac{u_{i+1}^n - u_i^n}{\Delta x} \end{align}$
(18)

Assuming $c > 0$ for this implementation, the upwind scheme update rule becomes:

\begin{align} u_i^{n+1} = u_i^n - \frac{c \Delta t}{\Delta x} \left( u_i^n - u_{i-1}^n \right) \end{align}

(19)

where:

$u_i^n$ is the numerical approximation of $u$ at spatial index $i$ and time level $n$ ,
$\Delta t$ and $\Delta x$ are the time and spatial step sizes, respectively.

The upwind scheme is actually simplier than the FTCS scheme. The following Python code implements it to solve the linear advection equation.

# Parameters

c  = 1.0    # advection speed
l  = 1.0    # domain size
dt = 0.005  # time step

nx = 100    # number of spatial points
nt = 200    # number of time steps

X, dx = np.linspace(0, l, nx, endpoint=False, retstep=True)  # spatial grid
U0    = np.sin(2*np.pi * X)  # initial condition

# Upwind scheme

def upwind(c, U0, dx, dt, n):
    sigma = c * dt / dx
    U = [U0]
    for _ in range(n):
        U0 = U[-1]
        U1 = U0 - sigma * (U0 - np.roll(U0,1))
        U.append(U1)
    return np.array(U)

U       = np.sin(2*np.pi * (X - c*dt*nt))  # analytical solution
Uupwind = upwind(c, U0, dx, dt, nt)        # numerical solution

plt.plot(X, U0,                label='Initial Condition')
plt.plot(X, U, ':',            label='Exact Solution')
plt.plot(X, Uupwind[-1], '.-', label='Upwind Scheme')
plt.xlabel('x')
plt.ylabel('u')
plt.legend()

anim = animate(X, Uupwind)

HTML(anim.to_html5_video())  # display animation
# anim.save('FTCS.mp4')        # save animation

Although the amplitude $u$ decrease as time evolve, the upwind scheme is at least stable!

# HANDSON: change `c`, `dt`, `dx`, or other parameters to observe how
#          the numerical solution changes.

# HANDSON: change the initial condition to, e.g., a boxcar function,
#          and observe how the numerical solution changes

Von Neumann Stability Analysis¶

Assuming again the numerical solution can be expressed as $u_i^n = G^n e^{ikx_i}$ , the upwind update equation becomes

\begin{align} G^{n+1} e^{ikx_i} = G^n e^{ikx_i} - \sigma \left( G^n e^{ikx_i} - G^n e^{ikx_{i-1}} \right) \end{align}

(20)

where $\sigma \equiv c \Delta t/\Delta x$ is again the Courant number.

Divide both sides by $G^n e^{ikx_i}$ ,

\begin{align} G = 1 - \sigma \left( 1 - e^{-ik\Delta x} \right). \end{align}

(21)

Using Euler’s formula and separate the real and imaginary parts

\begin{align} G = [1 - \sigma + \sigma \cos(k\Delta x)] - i \sigma \sin(k\Delta x). \end{align}

(22)

Calculate $|G|^2$ ,

\begin{align} |G|^2 &= [1 - \sigma + \sigma \cos(k\Delta x)]^2 + [\sigma \sin(k\Delta x)]^2 \\ &= (1 - \sigma)^2 + 2(1 - \sigma)\sigma\cos(k\Delta x) + \sigma^2\cos^2(k\Delta x) + \sigma^2 \sin^2(k\Delta x) \end{align}

(23)

Using $\cos^2\theta + \sin^2\theta = 1$ ,

\begin{align} |G|^2 = 1 - 2\sigma + 2\sigma^2 + 2(1 - \sigma)\sigma\cos(k\Delta x). \end{align}

(24)

We can easily plot it:

Kdx = np.linspace(-np.pi, np.pi, 65)

def absG(sigma):
    return np.sqrt(1 - 2 * sigma + 2 * sigma*sigma + 2*(1-sigma)*sigma*np.cos(Kdx))

for sigma in [0.5, 1, 1.5]:
    plt.plot(Kdx, absG(sigma), label=r"$\sigma={}$".format(sigma))

plt.xlabel(r"$k\Delta x$")
plt.ylabel(r"|G|")
plt.legend()

Clearly, the maximum of $|G|$ depends on $\sigma$ .

For $\sigma <= 1$ , consider $\cos(k\Delta x) = 1$ so that

\begin{align} |G|^2_{\text{max}} = 1. \end{align}

(25)

For $\sigma > 1$ , consider $\cos(k\Delta x) = -1$ so that

\begin{align} |G|^2_{\text{max}} = |1 - 2\sigma|^2 > 1. \end{align}

(26)

Therefore, the upwind scheme is stable as provided that the Courant number satisfies:

\begin{align} \sigma = \frac{c \Delta t}{\Delta x} \leq 1. \end{align}

(27)

# HANDSON: given what we know about stability now, change `c`, `dt`,
#          `dx`, or other parameters to observe how the numerical
#          solution changes.

Modified Equation Method¶

We observed the amplitude of the wave decreases when using the upwind scheme. This decay in amplitude is resemble of viscosity: in physical systems, viscosity acts to dissipate energy and smooth out oscillations. Similarly, many stable numerical schemes effectively remove “energy” from the discrete system in order to maintain stability. This is not surprising that the upwind method damps out $u(x,t)$ .

But how exactly does this damping arise from the discretization? To answer this, we turn to another powerful tool called the Modified Equation Method.

Whereas von Neumann stability analysis focuses on whether solutions grow or decay in time, modified equation analysis asks a different question:

What continuous differential equation is the numerical method really solving, up to leading-order truncation error?

By expanding the discrete update scheme in a Taylor series, we can write down an equivalent modified PDE. This PDE makes the numerical errors explicit, showing whether the scheme introduces artificial diffusion (numerical viscosity) or artificial dispersion (phase errors).

Modified Equation of Upwind Advection¶

Express the Upwind Scheme Using Taylor Series Expansions
Start with the upwind update equation:
$\begin{align} u_i^{n+1} = u_i^n - \sigma \left( u_i^n - u_{i-1}^n \right), \end{align}$
(28)
in addition to the grid point that we want to update, $u_i^n$ , there are two additional grid points. They are $u_i^{n+1}$ and $u_{i-1}^n$ .
Expand $u_i^{n+1}$ and $u_{i-1}^n$ around the point $(x_i, t^n)$ using Taylor series:
$\begin{align} u_i^{n+1} &= u(x_i, t^n) + \Delta t \left.\frac{\partial u}{\partial t}\right|_{x=x_i,t=t_n} + \frac{\Delta t^2}{2} \left.\frac{\partial^2 u}{\partial t^2}\right|_{x=x_i,t=t_n} + \mathcal{O}(\Delta t^3) \\ u_{i-1}^n &= u(x_i, t^n) - \Delta x \left.\frac{\partial u}{\partial x}\right|_{x=x_i,t=t_n} + \frac{\Delta x^2}{2} \left.\frac{\partial^2 u}{\partial x^2}\right|_{x=x_i,t=t_n} + \mathcal{O}(\Delta x^3) \end{align}$
(29)
Using the shorthands $\dot{u} \equiv \partial u/\partial t$ and $u' \equiv \partial u/\partial x$ , the above series can be written as
$\begin{align} u_i^{n+1} &= u_i^n + \Delta t\,\dot{u}_i^n + \frac{\Delta t^2}{2}\,\ddot{u}_i^n + \mathcal{O}(\Delta t^3) \\ u_{i-1}^n &= u_i^n - \Delta x\,{u' }_i^n + \frac{\Delta x^2}{2}\,{u'' }_i^n + \mathcal{O}(\Delta x^3) \end{align}$
(30)

Substitute the Taylor Expansions into the Upwind Scheme
Substitute the expansions into the upwind update equation:
$\begin{align} u_i^n + \Delta t\,\dot{u}_i^n + \frac{\Delta t^2}{2}\,\ddot{u}_i^n \approx u_i^n - \sigma \left[ u_i^n - \left( u_i^n - \Delta x\,{u'}_i^n + \frac{\Delta x^2}{2}\,{u''}_i^n \right) \right] \end{align}$
(31)
Because the derivatives are exact (instead of numerical), we can drop the functions’ evaluations at the dicrete points $x_i$ and $t^n$ . Simplify the equation, we have:
$\begin{align} \Delta t \frac{\partial u}{\partial t} + \frac{\Delta t^2}{2} \frac{\partial^2 u}{\partial t^2} \approx -\sigma \left( \Delta x \frac{\partial u}{\partial x} - \frac{\Delta x^2}{2} \frac{\partial^2 u}{\partial x^2} \right) \end{align}$
(32)

Rearrange and Substitute the Original PDE
Taking the spatial and temporal derivatives of the original advection equation $\partial u/\partial t = -c \partial u/\partial x$ , we have
$\begin{align} \frac{\partial^2 u}{\partial t^2} = -c \frac{\partial^2 u}{\partial x\partial t}, \quad \frac{\partial^2 u}{\partial t\partial x} = -c \frac{\partial^2 u}{\partial x^2} \end{align}$
(33)
Combine, we obtain the two-way wave equation:
$\begin{align} \frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}. \end{align}$
(34)
Substituting the wave equation into the modified equation, we obtain:
$\begin{align} \Delta t \frac{\partial u}{\partial t} + \frac{c^2\Delta t^2}{2} \frac{\partial^2 u}{\partial x^2} \approx -\sigma \left( \Delta x \frac{\partial u}{\partial x} - \frac{\Delta x^2}{2} \frac{\partial^2 u}{\partial x^2} \right). \end{align}$
(35)

Combine Like Terms
Put back the definition of $\sigma$ and rearrange,
$\begin{align} \frac{\partial u}{\partial t} + \frac{c^2\Delta t}{2} \frac{\partial^2 u}{\partial x^2} &\approx -c \left( \frac{\partial u}{\partial x} - \frac{\Delta x}{2} \frac{\partial^2 u}{\partial x^2} \right) \\ \frac{\partial u}{\partial t} +c \frac{\partial u}{\partial x} &\approx \frac{c}{2} (\Delta x - c\Delta t) \frac{\partial^2 u}{\partial x^2} \end{align}$
(36)
Define
$\begin{align} \nu_\text{upwind} \equiv c(\Delta x - c\Delta t) / 2 \end{align}$
(37)
be the numerical diffusion coefficient, the above equation reduces to
$\begin{align} \frac{\partial u}{\partial t} +c \frac{\partial u}{\partial x} \approx \nu_\text{upwind} \frac{\partial^2 u}{\partial x^2}. \end{align}$
(38)

From the modified equation analysis, we found that the Upwind scheme introduces an additional term of the form

\begin{align} \nu_{\text{upwind}} u_{xx}, \end{align}

(39)

where $\nu_{\text{upwind}} \propto (\Delta x - c \Delta t)$ depends on both the grid spacing and the chosen time step.

This term looks exactly like a viscous diffusion term, even though the original advection equation contains no viscosity at all. Hence the interpretation: the upwind scheme introduces “numerical viscosity/diffusion” to suppress instabilities.

Why numerical viscosity/diffusion matters?

Stabilizing Effect: Artificial viscosity is not necessarily a flaw. It is precisely what makes the Upwind Scheme stable, by damping out non-physical oscillations that can arise from discretization errors.
Trade-off with Accuracy: The downside is that this numerical viscosity smears out sharp features, such as wavefronts, contact discontinuities, or shocks. While this can improve stability, it comes at the cost of accuracy, especially in problems where sharp interfaces are physically meaningful.
Quantification via Modified Equation Analysis: The modified equation allows us to measure how much numerical diffusion is introduced and compare across schemes.

Thus, modified equation analysis connects stability theory (von Neumann analysis) and physical interpretation, showing how schemes balance between damping (diffusion) and oscillations (dispersion).

The artificial viscosity term also highlights the importance of the CFL condition:

\begin{align} \sigma = \frac{a \Delta t}{\Delta x} \leq 1. \end{align}

(40)

When $\sigma \leq 1$ , the artificial viscosity remains positive, ensuring stability.
When $\sigma > 1$ , the coefficient becomes negative, meaning the scheme adds anti-diffusion rather than diffusion, which destabilize the solution.

Therefore, the CFL condition directly controls the magnitude and sign of the artificial viscosity introduced by the scheme.

Modified Equation of FTCS Advection¶

To determine how the FTCS scheme modifies the original linear advection equation by introducing additional terms that represent numerical errors, specifically focusing on artificial diffusion or dispersion.

Start with the FTCS Update Equation
Recalling the FTCS scheme for the linear advection equation is given by:
$\begin{align} u_i^{n+1} = u_i^n - \frac{c \Delta t}{2 \Delta x} \left( u_{i+1}^n - u_{i-1}^n \right) \end{align}$
(41)
where:
- $u_i^n$ is the numerical approximation of $u$ at spatial index $i$ and time level $n$ ,
- $c$ is the constant advection speed,
- $\Delta t$ and $\Delta x$ are the time and spatial step sizes, respectively,
- $\sigma = c \Delta t/\Delta x$ is the Courant number.

Expand the Temporal and Spatial Terms Using Taylor Series
Expand $u_i^{n+1}$ around $(x_i, t^n)$ :
```
Undefined control sequence: \dddot at position 126: …Delta t^3}{6} {\̲d̲d̲d̲o̲t̲{u}\,}_i^n + \m…

\begin{align}
     u_i^{n+1} = u_i^n + \Delta t \dot{u}_i^n
     + \frac{\Delta t^2}{2} \ddot{u}_i^n + \frac{\Delta t^3}{6} {\dddot{u}\,}_i^n + \mathcal{O}(\Delta t^3)
   \end{align}
```
Expand $u_{i+1}^n$ and $u_{i-1}^n$ around $(x_i, t^n)$ :
$\begin{align} u_{i\pm1}^n &= u_i^n \pm \Delta x {u'}_i^n + \frac{\Delta x^2}{2} {u''}_i^n \pm \frac{\Delta x^3}{6} {u'''}_i^n + \mathcal{O}(\Delta x^4) \end{align}$
(43)

Substitute the Taylor Expansions into the FTCS Scheme

Substitute the expansions into the FTCS update equation:

Undefined control sequence: \dddot at position 110: …Delta t^3}{6} {\̲d̲d̲d̲o̲t̲{u}\,}_i^n \\
 …

\begin{align}
     u_i^n &+ \Delta t \dot{u}_i^n + \frac{\Delta t^2}{2} \ddot{u}_i^n + \frac{\Delta t^3}{6} {\dddot{u}\,}_i^n \\
     &\approx u_i^n - \frac{c \Delta t}{2 \Delta x} \left[
       \left( u_i^n + \Delta x {u'}_i^n + \frac{\Delta x^2}{2} {u''}_i^n + \frac{\Delta x^3}{6} {u'''}_i^n \right) -
       \left( u_i^n - \Delta x {u'}_i^n + \frac{\Delta x^2}{2} {u''}_i^n - \frac{\Delta x^3}{6} {u'''}_i^n \right)
   \right]
   \end{align}

Simplify,

Undefined control sequence: \dddot at position 89: …Delta t^2}{6} {\̲d̲d̲d̲o̲t̲{u}\,}_i^n
    …

\begin{align}
     \dot{u}_i^n + \frac{\Delta t}{2}\ddot{u}_i^n + \frac{\Delta t^2}{6} {\dddot{u}\,}_i^n
     &\approx - c\left[ {u'}_i^n + \frac{\Delta x^2}{6} {u'''}_i^n \right]
   \end{align}

Substitute and Rearrange
Similar to the upwind screen, we recall that the advection equation implies the wave equation $\partial^2 u/\partial t^2 = c^2 \partial^2 u/\partial x^2$ .
Taking an additional time derivative, we obtain:
$\begin{align} \frac{\partial^3 u}{\partial t^3} = c^2 \frac{\partial^3 u}{\partial x^2 \partial t} = - c^3\frac{\partial^3 u}{\partial x^3}. \end{align}$
(46)
Substitute this into the modified equation, we obtain
$\begin{align} \dot{u}_i^n + c^2 \frac{\Delta t}{2} {u''}_i^n - c^3 \frac{\Delta t^2}{6} {u'''}_i^n &\approx - c\left[ {u'}_i^n + \frac{\Delta x^2}{6} {u'''}_i^n \right]. \end{align}$
(47)
Rearrange and drop the indices $i$ and $n$ , the final modified equation is:
$\begin{align} \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} &\approx - \frac{c}{2} (c \Delta t ) \frac{\partial^2 u}{\partial x^2} + \frac{c}{6} (c^2 \Delta t^2 - \Delta x^2) \frac{\partial^3 u}{\partial x^3} \end{align}$
(48)

The modified equation of the FTCS advection scheme has a second-order “anti-diffusion” term and a third-order “dispersive” term. These terms affect the numerical solution in several ways:

Numerical Instability: Unlike the Upwind Scheme, which introduces artificial diffusion to enhance stability, the FTCS scheme has an anti-diffusion term that is unstable.
Introduction of Dispersive Errors: The (higher order) term $\partial^3 u/\partial x^3$ represents a dispersive error that causes different wave components to travel at slightly different speeds. This leads to phase errors where waveforms become distorted over time, deviating from the exact solution.
Impact on Solution Accuracy: The introduced dispersive term does not counteract the amplification of numerical errors. Instead, it modifies the original equation in a way that can exacerbate inaccuracies, especially for higher-frequency components of the solution. Over time, these errors accumulate, leading to significant deviations from the true solution, as evidenced by the von Neumann stability analysis.
Reinforcement of Von Neumann Stability Findings: The modified equation analysis complements the von Neumann stability analysis by providing a physical interpretation of why the FTCS scheme is unstable. The introduced dispersive errors contribute to the amplification of numerical oscillations, aligning with the conclusion that $|G| > 1$ for any $\sigma > 0$ .
Practical Considerations: Practitioners must recognize that the FTCS scheme not only fails to maintain stability but also introduces errors that distort the solution without providing any compensatory benefits.

Consequently, the FTCS scheme is unsuitable for solving advection-dominated PDEs where both stability and accuracy in capturing wave propagation are essential.

Lax-Wendroff Scheme¶

First-order schemes like the upwind method are simple to implement but often introduce significant numerical diffusion, which can blur important features of the solution. To address these limitations, higher-order schemes have been developed to provide more accurate and reliable results.

The Lax-Wendroff Scheme is a second-order accurate finite difference method designed to solve hyperbolic PDEs, such as the linear advection equation. Unlike first-order methods, the Lax-Wendroff scheme incorporates both temporal and spatial derivatives up to the second order. This allows it to capture wave propagation more accurately while minimizing numerical diffusion and dispersion.

The main advantage of the Lax-Wendroff Scheme is its ability to maintain higher accuracy without compromising stability. By extending the Taylor series expansion to include second-order terms, the scheme reduces the smearing effect seen in first-order methods. Additionally, it better preserves the shape and speed of waves, making it suitable for problems where precise wave behavior is crucial.

We begin by expanding $u(x, t + \Delta t)$ in time around $t$ using a Taylor series up to second order:

\begin{align} u(x, t + \Delta t) = u(x, t) + \Delta t \frac{\partial u}{\partial t} + \frac{(\Delta t)^2}{2} \frac{\partial^2 u}{\partial t^2} + \mathcal{O}(\Delta t^3) \end{align}

(49)

Substitute the advection equation and the wave equation into the Taylor series expansion:

\begin{align} u(x, t + \Delta t) = u(x, t) - c \Delta t \frac{\partial u}{\partial x} + \frac{c^2 (\Delta t)^2}{2} \frac{\partial^2 u}{\partial x^2} + \mathcal{O}((\Delta t)^3) \end{align}

(50)

and approximate the spatial derivatives using centered finite differences:

\begin{align} \frac{\partial u}{\partial x} &\approx \frac{u_{i+1}^n - u_{i-1}^n}{2 \Delta x} \\ \frac{\partial^2 u}{\partial x^2} &\approx \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{(\Delta x)^2}, \end{align}

(51)

we obtain:

\begin{align} u_i^{n+1} = u_i^n - c \Delta t \left(\frac{u_{i+1}^n - u_{i-1}^n}{2 \Delta x}\right) + \frac{c^2 \Delta t^2}{2} \left(\frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{\Delta x^2}\right) + \mathcal{O}(\Delta t^3) \end{align}

(52)

Simplify the equation to isolate $u_i^{n+1}$ , we obtain the Lax-Wendroff scheme for the linear advection equation

\begin{align} u_i^{n+1} = u_i^n - \frac{c \Delta t}{2 \Delta x} (u_{i+1}^n - u_{i-1}^n) + \frac{c^2 \Delta t^2}{2 \Delta x^2} (u_{i+1}^n - 2u_i^n + u_{i-1}^n) \end{align}

(53)

# Parameters

c  = 1.0    # advection speed
l  = 1.0    # domain size
dt = 0.005  # time step

nx = 101    # number of spatial points
nt = 200    # number of time steps

X, dx = np.linspace(0, l, nx, endpoint=False, retstep=True)  # spatial grid
U0    = np.sin(2*np.pi * X)  # initial condition

# Lax-Wendroff scheme

def LW(c, U0, dx, dt, n):
    sigma = c * dt / dx
    U = [U0]
    for _ in range(n):
        U0 = U[-1]
        U1 = U0 - sigma       * (np.roll(U0,-1)        - np.roll(U0,1)) / 2 \
                + sigma*sigma * (np.roll(U0,-1) - 2*U0 + np.roll(U0,1)) / 2
        U.append(U1)
    return np.array(U)

U   = np.sin(2*np.pi * (X - c*dt*nt)) # analytical solution
ULW = LW(c, U0, dx, dt, nt)           # numerical solution

plt.plot(X, U0,            label='Initial Condition')
plt.plot(X, U,       ':',  label='Exact Solution')
plt.plot(X, ULW[-1], '.-', label='Lax-Wendroff')
plt.xlabel('x')
plt.ylabel('u')
plt.legend()

anim = animate(X, ULW)

HTML(anim.to_html5_video())  # display animation
# anim.save('FTCS.mp4')        # save animation

Von Neumann Stability Analysis¶

Following previous derivation, it is straightforward to show

\begin{align} G = 1 + \sigma^2[\cos(k\Delta x) - 1] - i \sigma\sin(k\Delta x) \end{align}

(54)

Therefore,

\begin{align} |G|^2 = 1 + 4\sigma^2\sin^4\left(\frac{k\Delta x}{2}\right)(\sigma^2 - 1) \end{align}

(55)

Hence, $|G| \leq 1$ for all $k\Delta x$ iff $|\sigma| \leq 1$ . This lead to the same CFL condition

\begin{align} |c| \frac{\Delta t}{\Delta x} \leq 1. \end{align}

(56)

# HANDSON: Perform the von Neumann stability analysis to the
#          Lax-Wendroff Scheme.

Modified Equation¶

Also following previous derivation, we can also show the modified equation be

\begin{align} \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} &\approx \frac{c}{6} (c^2 \Delta t^2 - \Delta x^2) \frac{\partial^3 u}{\partial x^3} \end{align}

(57)

Note that this is the same as the FTCS modified equation but without the antidiffusion term.

Note also that when $\sigma = 1$ , the $\partial^3 u/\partial x^3$ term vanishes at leading order. When $\sigma \ne 1$ , the $\partial^3 u/\partial x^3$ term create phase (dispersion) error. For $|\sigma| < 1$ , there is also weak dissipation from higher-order terms.

# HANDSON: Perform modified equation analysis to the Lax-Wendroff
#          Scheme.