Topic: Finite Difference Methods and Conservation Viewpoint for Linear Advection
By the end of this notebook, you should be able to:
We consider the linear advection equation on a 1D domain \(x \in [0,L]\): \[u_t + c\,u_x = 0,\] where \(c\) is a constant wave speed.
We will study how different finite difference schemes approximate this simple behavior.
We discretize the domain \(x \in [0,L]\) into \(N\) grid points with spacing \(\Delta x = L/N\), and we advance in time using time step \(\Delta t\).
For simplicity, we will use periodic boundary conditions, so that \[u_0^n = u_N^n, \quad \text{or equivalently } u_{j+N}^n = u_j^n.\]
This allows us to focus on the behavior of the numerical schemes without dealing with boundaries.
import numpy as np
import matplotlib.pyplot as plt
# Basic grid setup
L = 1.0 # domain length
N = 200 # number of grid points
dx = L / N
x = np.linspace(0, L, N, endpoint=False)
c = 1.0 # advection speed
Tfinal = 0.5 # final time
# Courant number
lambda_target = 0.8
dt = lambda_target * dx / c
nt = int(Tfinal / dt)
print(f"dx = {dx:.4f}, dt = {dt:.4f}, number of time steps = {nt}")dx = 0.0050, dt = 0.0040, number of time steps = 125We will use a smooth Gaussian bump as an initial condition: \[u_0(x) = \exp\big(-100(x - 0.3)^2\big).\]
With periodic boundary conditions, this bump will wrap around the domain as it moves.
def initial_condition(x):
return np.exp(-100.0 * (x - 0.3)**2)
u0 = initial_condition(x)
plt.figure(figsize=(6, 3))
plt.plot(x, u0)
plt.xlabel("x")
plt.ylabel("u(x,0)")
plt.title("Initial Condition: Gaussian Bump")
plt.grid(True)
plt.show()
We start from the PDE \[u_t + c\,u_x = 0.\]
A simple explicit scheme is Forward-Time, Backward-Space (FTBS): \[ \frac{u_j^{n+1} - u_j^n}{\Delta t} + c \frac{u_j^n - u_{j-1}^n}{\Delta x} = 0. \]
Solving for \(u_j^{n+1}\): \[ u_j^{n+1} = u_j^n - \lambda\,(u_j^n - u_{j-1}^n), \quad \text{where } \lambda = \frac{c\,\Delta t}{\Delta x}. \]
This scheme is first order in both time and space, and for \(c>0\) it is stable when \(0 < \lambda \leq 1\).
We now implement FTBS with periodic boundaries.
def ftbs_advection(u0, c, dx, dt, nt):
'''
Forward-Time, Backward-Space scheme for u_t + c u_x = 0
with periodic boundary conditions.
'''
N = len(u0)
u = u0.copy()
lambda_c = c * dt / dx
history = [u.copy()]
for n in range(nt):
# periodic index j-1 -> (j-1) % N
u_new = u - lambda_c * (u - np.roll(u, 1))
u = u_new
history.append(u.copy())
return np.array(history)
history_ftbs = ftbs_advection(u0, c, dx, dt, nt)
plt.figure(figsize=(6, 3))
plt.plot(x, u0, label="t=0")
plt.plot(x, history_ftbs[nt//2], label=f"t={dt*(nt//2):.3f}")
plt.plot(x, history_ftbs[-1], label=f"t={dt*nt:.3f}")
plt.xlabel("x")
plt.ylabel("u")
plt.title("FTBS Scheme for Linear Advection")
plt.legend()
plt.grid(True)
plt.show()
Consider the Forward-Time, Central-Space (FTCS) scheme: \[ \frac{u_j^{n+1} - u_j^n}{\Delta t} + c\,\frac{u_{j+1}^n - u_{j-1}^n}{2\,\Delta x} = 0. \]
Solving for \(u_j^{n+1}\): \[ u_j^{n+1} = u_j^n - \frac{\lambda}{2}\,(u_{j+1}^n - u_{j-1}^n). \]
It turns out this scheme is unconditionally unstable for pure advection (no diffusion).
Hint: Use np.roll(u, -1) for \(u_{j+1}\) and np.roll(u, 1) for \(u_{j-1}\).
# TRY IT 1: Implement FTCS for linear advection
def ftcs_advection(u0, c, dx, dt, nt):
'''
TODO: Implement the Forward-Time, Central-Space scheme
for u_t + c u_x = 0 with periodic BCs.
'''
N = len(u0)
u = u0.copy()
lambda_c = c * dt / dx
history = [u.copy()]
for n in range(nt):
# TODO: add the line below with the FTCS update:
# u_new =
u_new = u.copy() # placeholder
history.append(u_new.copy())
u = u_new
return np.array(history)
# After you implement ftcs_advection, uncomment these lines to test:
# history_ftcs = ftcs_advection(u0, c, dx, dt, nt)
# plt.figure(figsize=(6, 3))
# plt.plot(x, u0, label="t=0")
# plt.plot(x, history_ftcs[nt//2], label=f"t={dt*(nt//2):.3f}")
# plt.plot(x, history_ftcs[-1], label=f"t={dt*nt:.3f}")
# plt.xlabel("x")
# plt.ylabel("u")
# plt.title("FTCS Scheme for Linear Advection (Expected to be Unstable)")
# plt.legend()
# plt.grid(True)
# plt.show()Observation:
This provides a nice contrast with FTBS, which is stable for \(0 < \lambda \le 1\) when \(c>0\).
For the linear advection equation, \[u_t + c\,u_x = 0,\] we can rewrite it in conservation form: \[u_t + \partial_x F(u) = 0, \quad F(u) = c\,u.\]
Integrating over a cell \([x_{j-1/2}, x_{j+1/2}]\): \[ \frac{d}{dt}\int_{x_{j-1/2}}^{x_{j+1/2}} u(x,t)\,dx = F\big(u(x_{j-1/2},t)\big) - F\big(u(x_{j+1/2},t)\big). \]
We will return to this idea in Week 14 when we formally introduce finite volume methods.
The Courant number is \[\lambda = \frac{c\,\Delta t}{\Delta x}.\]
For FTBS with \(c>0\), stability requires \(0 < \lambda \leq 1\).
# TRY IT 2: Explore the effect of lambda on FTBS stability
def run_ftbs_with_lambda(lambda_val):
dt_local = lambda_val * dx / c
nt_local = int(Tfinal / dt_local)
u0_local = initial_condition(x)
history_local = ftbs_advection(u0_local, c, dx, dt_local, nt_local)
return history_local, dt_local, nt_local
# TODO:
# 1. Call run_ftbs_with_lambda for lambda_val in [0.4, 0.8, 1.0, 1.2].
# 2. Plot the solutions at the final time on the same figure.
# 3. Observe which values of lambda lead to stability or instability,
# and which lead to more/less numerical diffusion.In this notebook we:
Next: We will build directly on the conservation form and derive finite volume methods, which are designed to enforce conservation at the discrete level by approximating fluxes across cell boundaries.