One PDE, Three Paths#

Put a PDE on a computer and the first cut is where the unknown sits. Same equation. Point value, basis coefficient, or cell average — three different methods.

One 1D advection-diffusion equation is enough to see the split.

$$\frac{\partial u}{\partial t} + a\frac{\partial u}{\partial x} = \nu\frac{\partial^2 u}{\partial x^2}$$

The dropdown switches FDM, FEM, and FVM, and the degrees of freedom move among nodes, basis coefficients, and cell centers.

FDM lines the points up, while FVM parks one value in the middle of each cell. FEM keeps the coefficient on the node and lets the basis run across the element.

---

1. Finite Difference Method (FDM)#

Core Idea#

“

Approximate derivatives directly as differences.

This is the most intuitive approach. It approximates derivatives by performing a Taylor expansion of the function values at grid points (nodes).

Mathematical Origin: Taylor Expansion#

Taylor expansion at point $x_i$:

$$u(x_i + \Delta x) = u(x_i) + \Delta x \frac{\partial u}{\partial x}\bigg|_i + \frac{\Delta x^2}{2}\frac{\partial^2 u}{\partial x^2}\bigg|_i + \cdots$$

From this, difference approximations are derived:

Forward Difference:

$$\frac{\partial u}{\partial x}\bigg|_i \approx \frac{u_{i+1} - u_i}{\Delta x} + O(\Delta x)$$

Central Difference:

$$\frac{\partial u}{\partial x}\bigg|_i \approx \frac{u_{i+1} - u_{i-1}}{2\Delta x} + O(\Delta x^2)$$

Second Derivative:

$$\frac{\partial^2 u}{\partial x^2}\bigg|_i \approx \frac{u_{i+1} - 2u_i + u_{i-1}}{\Delta x^2} + O(\Delta x^2)$$

Application to Advection-Diffusion#

Using central differences:

$$\frac{du_i}{dt} + a\frac{u_{i+1} - u_{i-1}}{2\Delta x} = \nu\frac{u_{i+1} - 2u_i + u_{i-1}}{\Delta x^2}$$

This results in a system of coupled ODEs for the unknowns $u_i$ at the grid points.

Pros and Cons#

Pros:

• Simple concept and easy implementation
• Highly efficient on structured grids
• Easy to achieve high-order accuracy (compact schemes, spectral-like schemes)
• Clean matrix structure on orthogonal grids (band matrix)

Cons:

• Difficult to apply to unstructured grids - this is a critical limitation
• Grid generation for complex geometries is challenging
• Does not automatically satisfy conservation laws

Typical Application Areas#

• DNS/LES (Turbulence simulations based on orthogonal grids)
• Weather/Ocean models (Structured grids)
• Seismic wave propagation simulation
• High-precision calculations based on compact schemes

---

2. Finite Element Method (FEM)#

Core Idea#

“

Approximate the solution as a linear combination of basis functions and minimize the weighted residual.

FEM does not solve the differential equation directly. Instead, it transforms it into an integral form called the weak formulation.

Mathematical Origin: Weak Form#

Multiply the original PDE (strong form) by a test function $w$ and integrate:

$$\int_\Omega w \left(\frac{\partial u}{\partial t} + a\frac{\partial u}{\partial x} - \nu\frac{\partial^2 u}{\partial x^2}\right) dx = 0$$

Applying integration by parts to the diffusion term:

$$\int_\Omega w \frac{\partial u}{\partial t}\,dx + \int_\Omega w\, a\frac{\partial u}{\partial x}\,dx + \int_\Omega \nu\frac{\partial w}{\partial x}\frac{\partial u}{\partial x}\,dx = \left[\nu w \frac{\partial u}{\partial x}\right]_\Gamma$$

This is the weak form. Note the key changes:

• The original second derivative was required, but after integration by parts, only the first derivative is needed.
• Requirements for solution continuity are relaxed ($C^1 \to C^0$).
• Boundary conditions are naturally included (Neumann BC = right-hand side).

Galerkin Approximation#

Express the solution as a linear combination of basis functions $\phi_j$:

$$u^h(x, t) = \sum_{j=1}^{N} U_j(t)\,\phi_j(x)$$

In the Galerkin method, test functions and basis functions are chosen from the same space ($w = \phi_i$):

$$\sum_j \left(\int_\Omega \phi_i \phi_j\,dx\right) \frac{dU_j}{dt} + \sum_j \left(\int_\Omega \phi_i\, a\, \phi_j'\,dx + \int_\Omega \nu\, \phi_i'\, \phi_j'\,dx\right) U_j = 0$$

In matrix form:

$$\mathbf{M}\frac{d\mathbf{U}}{dt} + \mathbf{K}\mathbf{U} = \mathbf{f}$$

• $\mathbf{M}$: Mass matrix
• $\mathbf{K}$: Stiffness matrix
• $\mathbf{f}$: Boundary/source vector

Choice of Basis Functions#

The most common choice is elements based on Lagrange polynomials:

Pros and Cons#

Pros:

• Natural for unstructured grids - handles complex shapes with triangular/tetrahedral meshes
• Mathematically rigorous error estimation (a priori / a posteriori) is possible
• Pairs well with adaptive mesh refinement
• Flexibility to choose between p-refinement (increasing polynomial degree) and h-refinement (grid refinement)

Cons:

• Does not satisfy conservation laws directly (Standard Galerkin)
• Unstable in advection-dominated problems (oscillation) - requires stabilization like SUPG or GLS
• Requires mass matrix inversion (or mass lumping)
• Higher implementation complexity compared to FDM/FVM

Typical Application Areas#

• Solid mechanics / Structural analysis (FEM's original home)
• Electromagnetic field analysis
• Heat transfer
• Biomechanics
• Geomechanics

---

3. Finite Volume Method (FVM)#

Core Idea#

“

Integrate conservation laws over control volumes and balance the fluxes passing through volume boundaries.

FVM reflects physical conservation laws most directly.

Mathematical Origin: Integral Conservation Law#

Write the advection-diffusion equation in integral form:

$$\frac{d}{dt}\int_{V_i} u\,dx + \oint_{\partial V_i} \left(au - \nu\frac{\partial u}{\partial x}\right) \cdot n\, dS = 0$$

For a control volume $V_i = [x_{i-1/2},\, x_{i+1/2}]$:

$$\frac{d\bar{u}_i}{dt} = -\frac{1}{\Delta x}\left(\hat{F}_{i+1/2} - \hat{F}_{i-1/2}\right)$$

Where $\bar{u}_i$ is the cell-averaged value and $\hat{F}_{i+1/2}$ is the numerical flux at the cell boundary.

Determining Numerical Fluxes#

The core of FVM is how to calculate cell boundary fluxes:

Advection term - based on Riemann solvers or upwind schemes:

$$\hat{F}^{adv}_{i+1/2} = \text{RiemannSolver}(\bar{u}_i,\, \bar{u}_{i+1})$$

Diffusion term - central difference:

$$\hat{F}^{diff}_{i+1/2} = -\nu\frac{\bar{u}_{i+1} - \bar{u}_i}{\Delta x}$$

Reconstruction#

Cell averages alone provide only first-order accuracy. For higher-order accuracy, the distribution inside the cell must be reconstructed:

MUSCL (2nd order):

$$u_{i+1/2}^L = \bar{u}_i + \frac{1}{2}\psi(r_i)\,(\bar{u}_i - \bar{u}_{i-1})$$

Where $\psi$ is a slope limiter and $r_i$ is the ratio of successive gradients.

WENO (5th order):

Achieves 5th-order accuracy by nonlinearly weighted averaging of three 2nd-order polynomial candidates:

$$u_{i+1/2}^L = \sum_{k=0}^{2} \omega_k\, q_k(x_{i+1/2})$$

Weights $\omega_k$ are determined based on smoothness indicators $\beta_k$, automatically giving higher weight to smooth stencils near discontinuities.

Face diffusion in the note above is $-\nu(\bar{u}_{i+1}-\bar{u}_i)/\Delta x$, an orthogonal spacing. That difference is the normal derivative only when the normal lies on the grid line. A skewed face does not get its direction from the neighbor averages alone. The non-orthogonal diffusion correction is what closes the gap. Drop it and a conservative update can still push diffusion off-axis.

The forward difference carried $O(\Delta x)$. Central and the second-derivative stencil carried $O(\Delta x^2)$. Forward keeps the Taylor term $\Delta x^2/2$. Central subtracts $u_{i+1}$ and $u_{i-1}$ and cancels it. Halve the spacing and the forward error drops by about half. The $O(\Delta x^2)$ pieces drop by about a quarter. Cell averages alone, with no reconstruction, stay first order. The extra order is bought on the way up to the face.

The $1/2$ in MUSCL is the factor that lifts a cell average linearly to the face. $\psi(r_i)=1$ keeps the whole slope. $\psi=0$ sends the face value back to $\bar{u}_i$, and the flux is first order again. WENO adds only three candidates, $k=0,1,2$. Each $q_k$ is a quadratic. $\omega_k$ cuts weight where $\beta_k$ is large. "Fifth order on a smooth interval" is the order of that weighted sum. When $\psi$ falls to 0 the fifth-order formula is still on the page, but that face has dropped to first order. A P1 triangle has 3 nodes, P2 has 6, P3 has 10. Raising the degree from 1 to 3 pushes the unknowns per element from 3 past 10. In the weak form, $\mathbf{M}$ holds $\int \phi_i \phi_j\,dx$ and $\mathbf{K}$ packs $\phi_i \phi_j'$ with $\nu\,\phi_i'\phi_j'$. Positive $\nu$ leaves the diffusion block nearly symmetric. The advective speed $a$ leaves the antisymmetric piece.

Pros and Cons#

Pros:

• Conservation laws are exactly satisfied at the discrete level - a decisive advantage in CFD
• Strong at capturing shock waves and discontinuities
• Applicable to unstructured grids
• Flux-based thinking directly aligns with physical intuition

Cons:

• Achieving high-order accuracy is more difficult than in FEM (reconstruction becomes complex)
• Can be less efficient than FEM for pure diffusion equations
• Implementation of high-order reconstruction on unstructured grids is tricky

Typical Application Areas#

• Across CFD (OpenFOAM, ANSYS Fluent, SU2)
• Compressible/Incompressible flow
• Multiphase flow, combustion, reacting flows
• Parts of dynamical cores in weather/climate models

---

Key Comparison#

Mathematical Starting Points#

"What" is Discretized?#

This is the most fundamental difference:

• FDM: Discretizes the differential operator ($\partial/\partial x$)
• FEM: Discretizes the solution space
• FVM: Discretizes the flux in integral equations

Location of Unknowns#

• FDM: Point values at grid points (nodes)
• FEM: Point values at nodes (coefficients of basis functions)
• FVM: Cell-averaged values

Comprehensive Comparison#

---

Modern Techniques: Blurring the Boundaries#

Recently, hybrid methods combining the advantages of all three techniques have been actively researched:

Discontinuous Galerkin (DG)#

Combines FEM's basis functions with FVM's flux concept. It approximates the solution with polynomials within each cell (FEM) and calculates fluxes at cell boundaries using Riemann solvers (FVM).

$$\int_{V_i} w \frac{\partial u^h}{\partial t}\,dx - \int_{V_i} \frac{\partial w}{\partial x} F(u^h)\,dx + \hat{F}_{i+1/2} w_{i+1/2}^- - \hat{F}_{i-1/2} w_{i-1/2}^+ = 0$$

It satisfies conservation + high-order accuracy + unstructured grid capabilities. However, computational costs are high, and limiters are required near shocks.

Spectral Difference / Flux Reconstruction#

Similar to DG but solves in differential form without integration to increase efficiency. It pairs very well with GPU acceleration and is gaining attention for next-generation CFD solvers.

Meshless Methods (SPH, etc.)#

An approach that eliminates the grid entirely, approximating based on particles. Strong for free-surface flows and large deformation problems, but accuracy/consistency issues exist.

---

If the integral has to survive, drop the pure difference stencil#

A box-domain DNS on an orthogonal mesh can live on the central $O(\Delta x^2)$ stencil. Three nodes build the second derivative, and the matrix stays banded. Bend the geometry onto an unstructured mesh and that stencil no longer has a neighbor to name. The weak form then swallows the shape. Standard Galerkin's $\mathbf{M}$ and $\mathbf{K}$ still do not promise that face fluxes cancel cell by cell.

When a shock passes and the integral has to come back as the same number, the update is FVM. The face flux $\hat{F}_{i+1/2}$ is a loss for the left cell and a gain for the right, so interior faces vanish in the global sum. High order and conservation in one scheme means a polynomial inside the cell and a Riemann flux on the face, as in a discontinuous Galerkin discretization. Leave the limiter off and that hybrid still rings in front of a shock.