Mathematical foundations of Wave Equation & CFL Stability.
From acoustic sound propagation and seismic earthquake waves to electromagnetic radiation and gravitational ripples in spacetime, dynamic non-dissipative transport is governed by hyperbolic partial differential equations. In 1928, Richard Courant, Kurt Friedrichs, and Hans Lewy published their landmark paper establishing the foundational stability principle of computational science: the Courant-Friedrichs-Lewy (CFL) condition. For explicit finite difference schemes, numerical information propagates at grid velocity Δx/Δt; if the physical wave speed c exceeds this numerical speed (Courant number C = c Δt / Δx > 1), physical domain characteristics escape the numerical domain of dependence, triggering instantaneous, violent exponential instability. In this computational laboratory, we implement the explicit leapfrog scheme for the 1D wave equation, simulate d'Alembert pulse splitting, investigate numerical dispersion, and verify the critical CFL stability parameter C = 1.000.
Viva voce questions & analytical theory
What is the physical interpretation of the Courant-Friedrichs-Lewy (CFL) condition?
The CFL condition states that the numerical domain of dependence of a finite difference scheme must contain the physical domain of dependence of the underlying hyperbolic PDE. If C = c Δt / Δx > 1, the physical characteristic lines travel faster than numerical information can propagate across the grid, rendering the numerical scheme fundamentally blind to physical wave signals and causing catastrophic instability.
How does d'Alembert's formula represent the exact solution to the 1D wave equation?
Jean le Rond d'Alembert proved in 1747 that the general solution to ∂²u/∂t² = c² ∂²u/∂x² is u(x, t) = f(x - ct) + g(x + ct), representing two non-dispersive wave packets propagating in opposite directions at constant velocity c without changing shape.
What is the amplification factor of the leapfrog wave scheme under von Neumann analysis?
The amplification factor satisfies the quadratic equation λ² - 2 (1 - 2 C² sin²(θ/2)) λ + 1 = 0. When C ≤ 1, the discriminant is non-positive, yielding complex conjugate roots with magnitude |λ| = 1 identically. This means the scheme is neutrally stable (non-dissipative) and conserves energy for C ≤ 1.
What happens when the Courant number equals exactly C = 1.0?
When C = c Δt / Δx = 1.0, the numerical scheme becomes mathematically exact: numerical grid characteristics coincide perfectly with physical characteristics, eliminating all numerical dispersion and transporting wave packets with zero distortion.
What is numerical dispersion in wave simulations?
When C < 1, the numerical phase velocity c_{num}(k) depends on the spatial wavenumber k: high-frequency short waves propagate slower than physical speed c. Over long integration times, sharp wave fronts disperse into parasitic trailing ripples.
How are non-reflecting boundary conditions (absorbing boundaries) formulated?
At open boundaries where waves should exit without reflecting back into the computational domain, one applies one-way wave equations (Engquist-Majda conditions) such as ∂u/∂t + c ∂u/∂x = 0 at the right boundary x = L, or surrounds the domain with a Perfectly Matched Layer (PML).
Why is the wave equation called 'hyperbolic'?
Classified by the discriminant B² - 4AC of general second-order PDEs A u_{xx} + B u_{xt} + C u_{tt} + ... = 0, the wave equation has A = -c², B = 0, C = 1, so B² - 4AC = 4c² > 0. A positive discriminant defines hyperbolic PDEs, which are characterized by finite propagation speeds along real characteristic curves.
Algorithm precautions and stability rules
- Always verify that Courant number C = c Δt / Δx ≤ 1.000 before initiating time integration to prevent catastrophic numerical overflow.
- Implement the initial velocity step using second-order Taylor expansion u¹ = u⁰ + (C²/2) δ_x² u⁰ to prevent degrading global accuracy to first order.
- Ensure boundary conditions at fixed ends u(0, t) = u(L, t) = 0 are enforced at every time step to prevent boundary energy injection.
- Verify that floating point precision remains 64-bit IEEE 754 during time-marching across thousands of wave periods.
- When simulating open domains, implement absorbing boundary conditions to prevent artificial wave packet reflection from grid boundaries.
- Monitor total discrete energy E^n = 0.5 ∑ [((u^{n+1} - u^{n-1})/(2Δt))² + c² ((u_{i+1}^n - u_i^n)/Δx)²] to detect instability onset.
Theoretical limitations
- Explicit leapfrog is strictly conditionally stable; in multi-dimensional domains, the CFL limit tightens to c Δt √(1/Δx² + 1/Δy²) ≤ 1.
- Numerical dispersion distorts high-frequency wave packets over long distances, requiring spectral methods or high-order compact stencils.
- The leapfrog scheme requires storage of three consecutive time levels, increasing memory footprint compared to two-level schemes.