PDE Numerical Methods

Discretizing partial differential equations on spatial grids: finite differences, finite elements, and spectral methods for the heat equation, wave equation, and Laplace's equation.

PDE Classification

The three canonical types of second-order PDEs, each with distinct physical behavior and numerical requirements:

Parabolic (Heat)

\(u_t = \alpha u_{xx}\)

Diffusion, smoothing. Information propagates at infinite speed.

Hyperbolic (Wave)

\(u_{tt} = c^2 u_{xx}\)

Waves, propagation. Finite speed, characteristics.

Elliptic (Laplace)

\(u_{xx} + u_{yy} = 0\)

Steady-state, boundary value problem. No time dependence.

1. Finite Difference Approximations

Replace derivatives with differences on a discrete grid. From Taylor expansion:

Forward: \(u'(x) \approx \frac{u(x+h) - u(x)}{h} + O(h)\)

Central: \(u'(x) \approx \frac{u(x+h) - u(x-h)}{2h} + O(h^2)\)

Second derivative: \(u''(x) \approx \frac{u(x+h) - 2u(x) + u(x-h)}{h^2} + O(h^2)\)

2. Heat Equation: Explicit Method

The 1D heat equation \(u_t = \alpha u_{xx}\) with the FTCS (Forward Time, Central Space) scheme:

\(u_j^{n+1} = u_j^n + r(u_{j+1}^n - 2u_j^n + u_{j-1}^n)\)

where \(r = \alpha \Delta t / (\Delta x)^2\). Stable only if \(r \leq 1/2\) (von Neumann analysis).

1D Heat Equation: Explicit FTCS

Python

script.py61 lines

import numpy as np

def heat_explicit(Nx, Nt, alpha, L, T, u0_func):
    """Solve 1D heat equation with explicit FTCS scheme."""
    dx = L / (Nx - 1)
    dt = T / Nt
    r = alpha * dt / dx**2

x = np.linspace(0, L, Nx)
    u = u0_func(x)

# Store snapshots
    snapshots = [(0, u.copy())]
    snapshot_times = [0.25, 0.5, 0.75, 1.0]
    snap_idx = 0

for n in range(Nt):
        u_new = u.copy()
        u_new[1:-1] = u[1:-1] + r * (u[2:] - 2*u[1:-1] + u[:-2])
        # Dirichlet BC: u(0,t) = u(L,t) = 0
        u_new[0] = 0
        u_new[-1] = 0
        u = u_new

t = (n + 1) * dt
        if snap_idx < len(snapshot_times) and t >= snapshot_times[snap_idx] * T:
            snapshots.append((t, u.copy()))
            snap_idx += 1

return x, snapshots, r

# Initial condition: sine wave
L, T, alpha = 1.0, 0.1, 1.0
u0 = lambda x: np.sin(np.pi * x)

# Stable case
Nx, Nt = 50, 5000
x, snaps, r = heat_explicit(Nx, Nt, alpha, L, T, u0)
print(f"1D Heat Equation: u_t = {alpha}*u_xx, L={L}")
print(f"Grid: {Nx} points, {Nt} time steps")
print(f"Stability parameter r = {r:.4f} (must be <= 0.5)")
print(f"Status: {'STABLE' if r <= 0.5 else 'UNSTABLE!'}")

print(f"\nSolution profile at selected times:")
print(f"{'x':>6}", end="")
for t, _ in snaps:
    print(f"  {'t='+f'{t:.3f}':>10}", end="")
print()
print("-" * (6 + 12 * len(snaps)))

sample_idx = np.linspace(0, Nx-1, 12, dtype=int)
for i in sample_idx:
    print(f"{x[i]:6.2f}", end="")
    for _, u_snap in snaps:
        print(f"  {u_snap[i]:10.6f}", end="")
    print()

# Exact solution: u(x,t) = sin(pi*x) * exp(-pi^2 * alpha * t)
u_exact = np.sin(np.pi * x) * np.exp(-np.pi**2 * alpha * T)
u_final = snaps[-1][1]
print(f"\nMax error vs exact at t={T}: {np.max(np.abs(u_final - u_exact)):.2e}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

3. Crank-Nicolson: The Best of Both Worlds

Average the explicit and implicit schemes to get second-order accuracy in both space and time, while remaining unconditionally stable:

\(u_j^{n+1} - \frac{r}{2}(u_{j+1}^{n+1} - 2u_j^{n+1} + u_{j-1}^{n+1}) = u_j^n + \frac{r}{2}(u_{j+1}^n - 2u_j^n + u_{j-1}^n)\)

Implicit: requires solving a tridiagonal system at each time step. \(O(h^2 + (\Delta t)^2)\) accuracy.

Crank-Nicolson for Heat Equation

Python

script.py89 lines

import numpy as np

def heat_crank_nicolson(Nx, Nt, alpha, L, T, u0_func):
    """Crank-Nicolson scheme for 1D heat equation."""
    dx = L / (Nx - 1)
    dt = T / Nt
    r = alpha * dt / dx**2

x = np.linspace(0, L, Nx)
    u = u0_func(x)

# Interior points only (BC handled separately)
    N = Nx - 2

# Build tridiagonal matrices
    # Left side: (I + r/2 * A) u^{n+1} = (I - r/2 * A) u^n
    diag_main = np.ones(N) * (1 + r)
    diag_off = np.ones(N - 1) * (-r / 2)

# Thomas algorithm for tridiagonal solve
    def tridiag_solve(a, b, c, d):
        """Solve tridiagonal system using Thomas algorithm."""
        n = len(d)
        c_ = np.zeros(n)
        d_ = np.zeros(n)
        c_[0] = c[0] / b[0]
        d_[0] = d[0] / b[0]
        for i in range(1, n):
            m = a[i-1] / (b[i] - a[i-1] * c_[i-1]) if i < n else 0
            denom = b[i] - a[i-1] * c_[i-1]
            if i < n - 1:
                c_[i] = c[i] / denom
            d_[i] = (d[i] - a[i-1] * d_[i-1]) / denom
        x = np.zeros(n)
        x[-1] = d_[-1]
        for i in range(n - 2, -1, -1):
            x[i] = d_[i] - c_[i] * x[i+1]
        return x

snapshots = [(0, u.copy())]
    for n in range(Nt):
        # RHS: (I - r/2 * A) u^n
        rhs = np.zeros(N)
        u_int = u[1:-1]
        rhs = u_int.copy()
        rhs[0] += (r/2) * u_int[0] if N > 0 else 0
        # Apply RHS operator
        rhs_full = np.zeros(N)
        for j in range(N):
            rhs_full[j] = (1 - r) * u_int[j]
            if j > 0:
                rhs_full[j] += (r/2) * u_int[j-1]
            if j < N - 1:
                rhs_full[j] += (r/2) * u_int[j+1]

# Solve tridiagonal system
        u_int_new = tridiag_solve(diag_off, diag_main, diag_off, rhs_full)
        u[1:-1] = u_int_new
        u[0] = 0
        u[-1] = 0

t = (n + 1) * dt
        if n == Nt // 4 or n == Nt // 2 or n == 3 * Nt // 4 or n == Nt - 1:
            snapshots.append((t, u.copy()))

return x, snapshots, r

L, T, alpha = 1.0, 0.1, 1.0
u0 = lambda x: np.sin(np.pi * x)

# Crank-Nicolson can use LARGE time steps (unconditionally stable)
Nx, Nt = 50, 50  # Only 50 time steps vs 5000 for explicit!
x, snaps, r = heat_crank_nicolson(Nx, Nt, alpha, L, T, u0)

print(f"Crank-Nicolson: 1D Heat Equation")
print(f"Grid: {Nx} points, only {Nt} time steps")
print(f"Stability parameter r = {r:.4f}")
print(f"Status: UNCONDITIONALLY STABLE (any r)\n")

# Compare with exact
u_exact = np.sin(np.pi * x) * np.exp(-np.pi**2 * alpha * T)
u_final = snaps[-1][1]
print(f"Max error vs exact at t={T}: {np.max(np.abs(u_final - u_exact)):.2e}")

# Show efficiency
print(f"\nEfficiency comparison for same grid ({Nx} points):")
print(f"  Explicit FTCS: needs r <= 0.5 => Nt >= {int(alpha * T / (0.5 * (L/(Nx-1))**2)) + 1}")
print(f"  Crank-Nicolson: used only Nt = {Nt}")
print(f"  Speedup: ~{int(alpha * T / (0.5 * (L/(Nx-1))**2)) // Nt}x fewer steps")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

4. The CFL Condition

The Courant-Friedrichs-Lewy condition is the fundamental stability criterion for explicit methods on hyperbolic (wave) equations:

\(\Delta t \leq C \frac{\Delta x}{|v_{\max}|}\)

The Courant number \(C = |v| \Delta t / \Delta x \leq 1\). Physical interpretation: the numerical domain of dependence must contain the physical domain of dependence.

Wave Equation: CFL Condition in Action

Python

script.py57 lines

import numpy as np

def wave_explicit(Nx, c, CFL, T, u0_func, v0_func=None):
    """Solve 1D wave equation u_tt = c^2 u_xx with explicit scheme."""
    L = 1.0
    dx = L / (Nx - 1)
    dt = CFL * dx / abs(c)
    Nt = int(T / dt) + 1
    r = (c * dt / dx)**2  # Should equal CFL^2

x = np.linspace(0, L, Nx)
    u_prev = u0_func(x)
    u_curr = u0_func(x)  # first step uses v0 = 0

# First time step (using initial velocity = 0)
    u_curr[1:-1] = u_prev[1:-1] + 0.5 * r * (u_prev[2:] - 2*u_prev[1:-1] + u_prev[:-2])

max_vals = [np.max(np.abs(u_curr))]
    stable = True

for n in range(2, Nt + 1):
        u_next = np.zeros(Nx)
        u_next[1:-1] = (2*u_curr[1:-1] - u_prev[1:-1] +
                         r * (u_curr[2:] - 2*u_curr[1:-1] + u_curr[:-2]))
        u_next[0] = 0
        u_next[-1] = 0

if np.max(np.abs(u_next)) > 1e10:
            stable = False
            break

max_vals.append(np.max(np.abs(u_next)))
        u_prev = u_curr
        u_curr = u_next

return stable, max_vals, Nt

c = 1.0
Nx = 100
u0 = lambda x: np.sin(2 * np.pi * x)

print("CFL Condition: 1D Wave Equation u_tt = c^2 u_xx")
print(f"Wave speed c = {c}, Grid points = {Nx}\n")

print(f"{'CFL':>6}  {'Stable?':>8}  {'Steps':>8}  {'Max |u| final':>15}")
print("-" * 42)
for CFL in [0.5, 0.8, 0.95, 1.0, 1.01, 1.1, 1.5]:
    stable, max_vals, Nt = wave_explicit(Nx, c, CFL, T=2.0, u0_func=u0)
    max_u = max_vals[-1] if stable else float('inf')
    status = "YES" if stable else "BLOWUP"
    max_str = f"{max_u:.6f}" if stable else "INFINITY"
    print(f"{CFL:6.2f}  {status:>8}  {Nt:8d}  {max_str:>15}")

print(f"\nCFL <= 1.0: STABLE (max|u| stays bounded)")
print(f"CFL >  1.0: UNSTABLE (solution blows up)")
print(f"\nPhysical meaning: numerical info must travel")
print(f"at least as fast as physical waves.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

5. Finite Element Method: Basics

Instead of approximating derivatives (finite differences), FEM approximates the solution itself as a sum of basis functions\(u(x) \approx \sum_j c_j \phi_j(x)\) and finds coefficients via the weak (variational) form.

Weak Form of Poisson's Equation

Strong form: \(-u''(x) = f(x)\) on \([0, 1]\), \(u(0) = u(1) = 0\)

Weak form: Find \(u \in H_0^1\) such that for all test functions \(v \in H_0^1\):

\(\int_0^1 u'(x) v'(x) \, dx = \int_0^1 f(x) v(x) \, dx\)

Integration by parts trades one derivative from \(u\) to \(v\), weakening the smoothness requirement. With piecewise-linear (hat function) basis, this yields a tridiagonal stiffness matrix.

1D Finite Element Method for Poisson Equation

Python

script.py57 lines

import numpy as np

def fem_poisson_1d(N, f_func):
    """Solve -u'' = f on [0,1] with u(0)=u(1)=0 using linear FEM."""
    h = 1.0 / N
    x = np.linspace(0, 1, N + 1)

# Stiffness matrix K and load vector F
    # For linear hat functions:
    # K_ij = integral(phi_i' * phi_j') = 1/h on diagonal, -1/h off-diagonal
    # F_i = integral(f * phi_i) ~ h * f(x_i)  (midpoint rule)

# Interior nodes: 1 to N-1
    n = N - 1
    K = np.zeros((n, n))
    F = np.zeros(n)

for i in range(n):
        K[i, i] = 2.0 / h
        if i > 0:
            K[i, i-1] = -1.0 / h
        if i < n - 1:
            K[i, i+1] = -1.0 / h
        F[i] = h * f_func(x[i + 1])

# Solve
    u_int = np.linalg.solve(K, F)

# Full solution with BCs
    u = np.zeros(N + 1)
    u[1:-1] = u_int
    return x, u

# Test: -u'' = pi^2 * sin(pi*x), exact: u = sin(pi*x)
f = lambda x: np.pi**2 * np.sin(np.pi * x)
u_exact_func = lambda x: np.sin(np.pi * x)

print("1D FEM: -u'' = pi^2 sin(pi x), u(0)=u(1)=0")
print(f"Exact solution: u(x) = sin(pi x)\n")

print(f"{'N':>6}  {'h':>10}  {'Max error':>12}  {'Order':>6}")
print("-" * 40)
prev_err = None
for N in [10, 20, 40, 80, 160]:
    x, u = fem_poisson_1d(N, f)
    u_exact = u_exact_func(x)
    err = np.max(np.abs(u - u_exact))
    if prev_err is not None:
        order = np.log2(prev_err / err)
        print(f"{N:6d}  {1.0/N:10.4f}  {err:12.4e}  {order:6.2f}")
    else:
        print(f"{N:6d}  {1.0/N:10.4f}  {err:12.4e}  {'---':>6}")
    prev_err = err

print("\nOrder ~2.0: linear FEM gives O(h^2) convergence")
print("(Same as finite differences for this problem, but FEM")
print(" generalizes to complex geometries and higher dimensions)")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

6. Spectral Methods

Use global basis functions (Fourier modes, Chebyshev polynomials) instead of local ones. For smooth solutions, spectral methods converge exponentially (faster than any polynomial rate).

Spectral Method: Solving Poisson with FFT

Python

script.py44 lines

import numpy as np

def spectral_poisson(N, f_func):
    """Solve -u'' = f on [0, 2*pi] periodic using FFT."""
    x = np.linspace(0, 2*np.pi, N, endpoint=False)
    f_vals = f_func(x)

# FFT of RHS
    f_hat = np.fft.fft(f_vals)

# In Fourier space: -u'' => k^2 * u_hat = f_hat
    k = np.fft.fftfreq(N, d=1.0/N)
    k2 = k**2

# Solve in Fourier space (skip k=0 mode)
    u_hat = np.zeros(N, dtype=complex)
    u_hat[1:] = f_hat[1:] / k2[1:]
    u_hat[0] = 0  # set mean to zero

# Inverse FFT
    u = np.real(np.fft.ifft(u_hat))
    return x, u

# Test: -u'' = sin(x) on [0, 2pi] periodic
# Exact: u = sin(x) (since sin'' = -sin)
f = lambda x: np.sin(x)
u_exact_func = lambda x: np.sin(x)

print("Spectral Method: -u'' = sin(x), periodic BC")
print(f"Exact: u(x) = sin(x)\n")

print(f"{'N':>6}  {'Max error':>14}  {'Notes':>20}")
print("-" * 44)
for N in [4, 8, 16, 32, 64, 128]:
    x, u = spectral_poisson(N, f)
    u_exact = u_exact_func(x)
    err = np.max(np.abs(u - u_exact))
    note = "machine precision" if err < 1e-13 else ""
    print(f"{N:6d}  {err:14.2e}  {note:>20}")

print("\nSpectral methods reach machine precision with very few points!")
print("For smooth periodic problems, exponential convergence")
print("far outperforms algebraic convergence (O(h^p)).")
print("\nCost: O(N log N) per solve using FFT")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

7. Laplace Equation: 2D Steady State

2D Laplace Equation with Jacobi Iteration

Python

script.py48 lines

import numpy as np

def laplace_2d_jacobi(Nx, Ny, bc_func, tol=1e-6, max_iter=10000):
    """Solve Laplace equation on unit square using Jacobi iteration."""
    dx = 1.0 / (Nx - 1)
    dy = 1.0 / (Ny - 1)
    u = np.zeros((Ny, Nx))

# Set boundary conditions
    x = np.linspace(0, 1, Nx)
    y = np.linspace(0, 1, Ny)
    u[-1, :] = bc_func(x)  # top boundary: u(x, 1)
    # Other boundaries = 0

for iteration in range(max_iter):
        u_old = u.copy()
        u[1:-1, 1:-1] = 0.25 * (u_old[2:, 1:-1] + u_old[:-2, 1:-1] +
                                   u_old[1:-1, 2:] + u_old[1:-1, :-2])
        # Check convergence
        diff = np.max(np.abs(u - u_old))
        if diff < tol:
            return u, x, y, iteration + 1

return u, x, y, max_iter

# BC: u(x,1) = sin(pi*x), all other boundaries = 0
# Exact: u(x,y) = sin(pi*x) * sinh(pi*y) / sinh(pi)
bc = lambda x: np.sin(np.pi * x)

N = 30
u, x, y, iters = laplace_2d_jacobi(N, N, bc)

print(f"2D Laplace Equation: u_xx + u_yy = 0")
print(f"Grid: {N}x{N}, converged in {iters} Jacobi iterations\n")

# Compare with exact solution
X, Y = np.meshgrid(x, y)
u_exact = np.sin(np.pi * X) * np.sinh(np.pi * Y) / np.sinh(np.pi)
max_err = np.max(np.abs(u - u_exact))
print(f"Max error vs exact: {max_err:.4e}")

# Show solution along middle row
print(f"\nSolution at y = 0.5:")
mid_j = N // 2
print(f"{'x':>6}  {'u_numeric':>12}  {'u_exact':>12}  {'|error|':>12}")
print("-" * 46)
for i in range(0, N, max(1, N//10)):
    print(f"{x[i]:6.2f}  {u[mid_j,i]:12.6f}  {u_exact[mid_j,i]:12.6f}  {abs(u[mid_j,i]-u_exact[mid_j,i]):12.2e}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server