Partial Differential Equations: FDM, FEM, and the CFL Condition

Von Neumann, Friedrichs, and the Birth of Numerical Stability

John von Neumann and Kurt Friedrichs were running numerical simulations of shock waves for US defense calculations. The equations were correct. The grid was fine. But the results exploded - numbers diverged to infinity. The problem was not the physics and not the code: the time step was too large relative to the spatial step. Their 1950 stability analysis produced a mathematical condition that now bears their names alongside Courant: the CFL criterion. Every climate model, aerodynamics solver, plasma simulator, and seismic code checks this condition before each run.

Finite Difference Methods and the CFL Condition

1950. Von Neumann and Friedrichs are running numerical simulations of shock waves for US defense calculations. The equations are correct. The grid is reasonable. But the solution blows up - numbers spiral to infinity after a few iterations. The problem is not the physics and not the code: the time step is too large relative to the spatial step. Their stability analysis gave a mathematical condition now known as the CFL criterion. Every climate model, aerodynamics solver, and plasma simulator on the planet checks this condition before running.

A PDE - partial differential equation - describes a field: temperature at every point of a plate, pressure at every point in the atmosphere, probability density across phase space. The general form: $\frac{\partial u}{\partial t} = \mathcal{L}[u]$, where $\mathcal{L}$ is a spatial operator (Laplacian, gradient, divergence). Analytical solutions exist for a handful of special cases. Numerical solutions are needed everywhere else.

The Finite Difference Method (FDM) is the most direct idea: replace derivatives with difference quotients. Space $x \in [0, L]$ is divided into $N$ nodes with step $\Delta x = L/N$. Time is discretized into steps $\Delta t$. Three approximations for the first derivative exist.

For the heat equation $\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}$, the explicit FTCS scheme (Forward Time, Central Space) uses current-layer values to compute the next layer. Simple. Cheap. But conditionally stable.

**CFL Condition (Courant-Friedrichs-Lewy):** for the explicit heat equation scheme: $r = \frac{\alpha \Delta t}{\Delta x^2} \leq \frac{1}{2}$. For the wave equation: $c \frac{\Delta t}{\Delta x} \leq 1$. Violation causes errors to grow exponentially. This is not a warning. This is the exact boundary between a simulator and a garbage generator.

The Crank-Nicolson implicit scheme averages explicit and implicit layers. Result: unconditional stability - any $\Delta t$ works. But each time step requires solving a linear system. This is the fundamental trade-off: explicit schemes are cheap per step but time-step limited; implicit schemes are stable at any step size but more expensive per step. ECMWF global weather forecasts, oceanic circulation models - all implicit, because the time steps needed are enormous.

Finite Element Method (FEM)

FDM works perfectly on rectangular grids. But an aircraft wing is not a rectangle. A human skull is not a cube. A turbine blade is not a cylinder. The Finite Element Method (FEM) was designed precisely for irregular geometry - and became the industry standard for structural and fluid analysis.

FEM takes a different route. Instead of approximating derivatives with differences, it uses a **variational formulation**. The PDE is multiplied by a test function $v$ and integrated over the domain: $\int_\Omega \frac{\partial u}{\partial t} v \, d\Omega = \int_\Omega \mathcal{L}[u] \, v \, d\Omega$. This is the weak (variational) form. The solution $u$ is approximated as a linear combination of basis functions defined on each element.

The domain is split into triangles (2D) or tetrahedra (3D). On each element, $u$ is approximated by a low-degree polynomial. The resulting system for the coefficients is a **matrix problem**: $\mathbf{K} \mathbf{u} = \mathbf{f}$, where $\mathbf{K}$ is the global stiffness matrix - sparse. This is what Ansys, Abaqus, OpenFOAM, and FEniCS do under the hood.

**FEM vs FDM:** FDM is simpler to implement on structured grids and faster for simple geometries. FEM handles complex geometry but has higher implementation overhead. In industry (aerodynamics, structural analysis, biomechanics) - FEM. In meteorology, oceanography, and fundamental physics - often FDM or FVM.

The stiffness matrix $\mathbf{K}$ is sparse: each node couples only to its neighbors. For a 3D mesh with $N$ nodes - $O(N)$ nonzeros out of $O(N^2)$ possible entries. This is what makes FEM scalable - problems with millions of degrees of freedom are solved with iterative methods (CG, GMRES with preconditioner) in near-linear time.

Explicit vs Implicit Schemes: Stability and Cost

Every climate simulator faces the same choice: explicit schemes give simple code and cheap steps, but demand small $\Delta t$ or CFL blows up. Implicit schemes are stable at any $\Delta t$, but each step requires solving a linear system. This is not an academic question - it determines how many millions of CPU-hours a forecast costs.

**Von Neumann stability analysis** is the universal tool for checking a scheme. The solution is expanded in a Fourier series: $u_j^n = \xi^n e^{ikj\Delta x}$. The scheme is stable if the amplification factor satisfies $|\xi| \leq 1$ for all wave numbers $k$. Here is what this looks like for the heat equation.

Operator splitting is an elegant compromise. A time step is split into stages: explicit for non-stiff terms (advection), then implicit for stiff terms (diffusion). This is IMEX (IMplicit-EXplicit) integration. WENO schemes for shock waves, spectral methods for DNS turbulence - all variations on this theme.

**Numbers at scale:** ECMWF runs its IFS model at ~9 km resolution, $\Delta t = 450$ s, 137 vertical levels. One 10-day forecast takes ~15 minutes on 1000+ cores. An explicit scheme with CFL constraint would need $\Delta t \approx 30$ s - 15x more steps, 15x the compute time.

In machine learning, the CFL analog appears in Physics-Informed Neural Networks (PINNs). The network minimizes the PDE residual at collocation points. No grid, no CFL - but a different instability: balancing residual losses between the equation and boundary conditions. DeepMind and NVIDIA Neural Operator (FNO, GNO) go further - training operators that map one field to another, bypassing step-by-step integration entirely.

Key Ideas

• **FDM** replaces derivatives with difference quotients on regular grids - simple, fast, limited to rectangular domains
• **CFL condition:** explicit schemes are stable only when $r = \alpha \Delta t / \Delta x^2 \leq 1/2$ for diffusion and $c \cdot \Delta t / \Delta x \leq 1$ for waves - violation gives exponential error growth
• **FEM** solves the variational form of the PDE on unstructured meshes - the industry standard for arbitrary geometries, reduces to a sparse linear system $\mathbf{K}\mathbf{u} = \mathbf{f}$
• **Explicit schemes** are cheap per step but time-step limited; **implicit schemes** are unconditionally stable but solve a linear system each step - real simulators use IMEX operator splitting
• **PINNs and neural operators** (FNO, GNO) bypass step-by-step integration by learning the solution operator directly from data

Related Topics

PDE methods bring together linear algebra, numerical analysis, and high-performance computing:

• ODEs and Runge-Kutta — ODEs are PDEs without spatial derivatives; the same stability questions apply
• Systems of Linear Equations — FEM reduces the PDE to a sparse SLE Ku=f solved by iterative methods
• High-Performance Computing — Large FDM/FEM problems require MPI domain decomposition and GPU acceleration
• Computational Physics — Direct application of FDM/FEM to fluid dynamics, electromagnetism, and elasticity

Вопросы для размышления

• The CFL condition says information cannot propagate faster than the numerical scheme allows. What is the physical interpretation behind this constraint?
• Why does halving $\Delta x$ in an explicit diffusion scheme require reducing $\Delta t$ by a factor of four - and how does this compound in 3D simulations?
• Physics-Informed Neural Networks avoid the CFL condition. What makes this possible, and what new stability problems arise in their place?

Связанные уроки

• sci-05 — ODEs are a special case of PDEs; Runge-Kutta stability is a prerequisite
• sci-03 — FEM reduces the PDE to a sparse linear system - SLE background needed
• sci-07 — Functional optimization is the next level after variational FEM formulations
• sci-08 — Stochastic PDEs and Monte Carlo methods solve problems beyond FDM reach
• sci-11 — Computational physics applies FDM/FEM directly to fluids, elasticity, EM fields
• nm-07