One PDE, Three Paths#

To solve partial differential equations (PDEs) on a computer, continuous space must be discretized. Even for the same equation, the resulting numerical technique becomes completely different depending on the philosophy used for discretization.

Let's use a simple 1D advection-diffusion equation as an example:

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

We will examine how FDM, FEM, and FVM each approach this single equation.

---

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.

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.

---

Conclusion: Which One Should You Choose?#

There is no "magic" technique. The nature of the problem dictates the choice:

• Structural analysis, heat transfer, electromagnetics: FEM
• Compressible flow, shocks, multiphase flow: FVM
• DNS, high-precision structured grid calculations: FDM
• Simultaneous high-order accuracy + conservation: DG (FEM + FVM hybrid)

The key is to understand the mathematical starting points and limitations of each technique. This allows you to choose the right tool for the problem and interpret the results correctly.

Switch between FDM/FEM/FVM in the dropdown to see how each places its degrees of freedom on the same grid.