PDE Solvers

Partial differential equations model phenomena where quantities depend on multiple independent variables (space, time). Numerical PDE methods are the backbone of scientific computing, simulation, and engineering.

Classification of PDEs

Second-Order Linear PDEs

General form in 2D: Au_xx + 2Bu_xy + Cu_yy + Du_x + Eu_y + Fu = G.

Classification by discriminant B² - AC:

| Type | Condition | Prototype | Physical Model | |---|---|---|---| | Elliptic | B² - AC < 0 | Laplace: u_xx + u_yy = 0 | Steady-state (equilibrium) | | Parabolic | B² - AC = 0 | Heat: u_t = α·u_xx | Diffusion, smoothing | | Hyperbolic | B² - AC > 0 | Wave: u_tt = c²·u_xx | Wave propagation |

Each type requires different numerical approaches and boundary conditions.

Boundary Conditions

• Dirichlet: u = g on boundary (value specified)
• Neumann: ∂u/∂n = g on boundary (derivative specified)
• Robin (mixed): αu + β∂u/∂n = g
• Periodic: u(0) = u(L)

Finite Difference Methods

Replace derivatives with finite differences on a grid.

Heat Equation: u_t = α·u_xx

Explicit (FTCS — Forward Time, Central Space):

u_i^{n+1} = u_i^n + r(u_{i+1}^n - 2u_i^n + u_{i-1}^n)

where r = αΔt/Δx².

Stability: Stable iff r ≤ 1/2 (CFL condition). This constrains Δt ≤ Δx²/(2α) — very restrictive for small Δx!

Implicit (Backward Euler):

u_i^{n+1} - r(u_{i+1}^{n+1} - 2u_i^{n+1} + u_{i-1}^{n+1}) = u_i^n

Requires solving a tridiagonal system at each time step. Unconditionally stable (any Δt).

Crank-Nicolson (average of explicit and implicit):

u_i^{n+1} = u_i^n + r/2(δ²u^{n+1} + δ²u^n)

where δ²u = u_{i+1} - 2u_i + u_{i-1}.

Unconditionally stable and second-order in both time and space: O(Δt² + Δx²). The standard choice for parabolic equations.

Wave Equation: u_tt = c²·u_xx

Central differences in both time and space:

u_i^{n+1} = 2u_i^n - u_i^{n-1} + r²(u_{i+1}^n - 2u_i^n + u_{i-1}^n)

where r = cΔt/Δx (Courant number).

CFL condition: r ≤ 1. The wave speed must not exceed the grid speed Δx/Δt.

Laplace/Poisson Equation: u_xx + u_yy = f

Five-point stencil:

(u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} - 4u_{i,j}) / h² = f_{i,j}

Results in a large sparse linear system Au = f. Solved by iterative methods (Gauss-Seidel, SOR, multigrid, conjugate gradient).

For an N × N grid: N² unknowns, O(N²) non-zeros in A.

Convergence Analysis

Consistency: The FD scheme approximates the PDE as Δx, Δt → 0.

Stability: Errors don't grow unboundedly. Analyzed by von Neumann analysis (Fourier modes).

Lax Equivalence Theorem: For a consistent scheme, stability ⟺ convergence. This is the fundamental theorem of numerical PDEs.

Finite Element Methods (FEM)

The most versatile method for complex geometries and unstructured meshes.

Weak Formulation

Multiply the PDE by a test function v and integrate by parts:

For -∇²u = f with Dirichlet BCs:

∫_Ω ∇u · ∇v dΩ = ∫_Ω fv dΩ    for all test functions v

This is the weak form (variational formulation). Requires less smoothness than the original PDE.

Galerkin Method

Approximate u ≈ Σⱼ cⱼ φⱼ(x) where {φⱼ} are basis functions on a mesh.

Substituting into the weak form gives a linear system: Kc = f.

Stiffness matrix: K_ij = ∫ ∇φᵢ · ∇φⱼ dΩ Load vector: f_i = ∫ f·φᵢ dΩ

K is sparse, symmetric positive definite (for elliptic problems).

Basis Functions

Linear elements (P1): Piecewise linear on triangles. "Hat functions." Each basis function is 1 at one node, 0 at all others.

Quadratic elements (P2): Piecewise quadratic. 6 nodes per triangle (vertices + edge midpoints). Better accuracy.

Higher-order: P3, P4, etc. Or use different element shapes (quads, hexahedra, tetrahedra).

Mesh Generation

• Triangulation: Delaunay triangulation, advancing front, octree-based
• Quality: Avoid thin/degenerate elements (poor condition number)
• Adaptive refinement: Refine elements where error is large (a posteriori error estimation)

Error Estimates

For P1 elements on a quasi-uniform mesh with element size h:

‖u - u_h‖_{L²} = O(h²)
‖u - u_h‖_{H¹} = O(h)

Higher-order elements give higher-order convergence.

Finite Volume Methods

Integrate the conservation law over each control volume:

∂u/∂t + ∇ · F(u) = 0

Integrate over cell Ωᵢ:

d/dt ∫_Ωᵢ u dΩ = -∮_{∂Ωᵢ} F · n dS

The flux through each face is approximated and balanced.

Advantages: Naturally conserves quantities. Handles discontinuities (shocks) well. Dominant in CFD (computational fluid dynamics).

Godunov-type methods: Solve Riemann problems at cell interfaces. MUSCL reconstruction for higher order.

Spectral Methods

Represent the solution as a sum of global basis functions (Fourier modes, Chebyshev polynomials).

u(x) = Σₖ ûₖ φₖ(x)

Fourier spectral: For periodic problems. Uses FFT. Converges exponentially fast for smooth solutions.

Chebyshev spectral: For non-periodic problems on [-1, 1]. Exponential convergence for smooth solutions.

Advantages: Spectral (exponential) convergence for smooth problems. Very high accuracy with few degrees of freedom.

Disadvantages: Poor for discontinuities (Gibbs phenomenon). Complex geometries require spectral element methods.

Multigrid Methods

The most efficient solver for elliptic and parabolic PDEs.

Problem

Iterative methods (Jacobi, Gauss-Seidel) quickly reduce high-frequency error but stall on low-frequency error.

Key Idea

Smooth error on the fine grid appears oscillatory on a coarser grid, where iterative methods are effective. Cycle between grids:

1. Relax (smooth) on fine grid: reduce high-frequency error
2. Restrict residual to coarse grid
3. Solve on coarse grid (recursively or exactly)
4. Prolongate (interpolate) correction back to fine grid
5. Relax again on fine grid

V-Cycle

Smooth → Restrict → [Coarse solve] → Prolongate → Smooth

The coarse solve is itself a V-cycle → recursion down to the coarsest grid.

Complexity

O(N) operations to solve an N-unknown elliptic problem! This is optimal — you can't even read the answer faster than O(N).

Compare: Direct solvers O(N^(3/2)) in 2D, iterative (CG) O(N^(3/2)) to O(N log N).

Algebraic Multigrid (AMG)

Automatically constructs coarse grids from the matrix structure (no geometric mesh needed). Used as a preconditioner for CG/GMRES.

Domain Decomposition

Split the domain into subdomains, solve on each, and couple at interfaces.

Schwarz methods: Overlapping subdomains. Alternate between solving on each subdomain with boundary data from the other.

Schur complement: Reduce to an interface problem. Solve the interface system, then recover the interior.

Parallel implementation: Each subdomain on a different processor. Communication only at interfaces. Scales to millions of cores.

Applications in CS

• Computational fluid dynamics: Weather prediction, aerodynamics (aircraft design), blood flow simulation. Navier-Stokes equations solved by FVM or FEM.
• Computer graphics: Fluid simulation (smoke, water, fire) for games and movies. Poisson equation for pressure. Heat equation for diffusion effects.
• Structural analysis: Stress, deformation, vibration of bridges, buildings, car bodies. FEM is the standard.
• Electromagnetic simulation: Antenna design, signal integrity, EMI. Maxwell's equations solved by FDTD or FEM.
• Machine learning: Physics-informed neural networks (PINNs) learn PDE solutions. Neural operators (FNO, DeepONet) approximate PDE solution operators.
• Image processing: Poisson image editing, inpainting, denoising (diffusion equation), level sets for segmentation.
• Semiconductor device simulation: Drift-diffusion equations, Poisson equation for electrostatics.
• Geophysics: Seismic wave propagation, reservoir simulation, mantle convection.