01 / NUMERICAL VERIFICATION

Can We Trust FEM?

PRECOMPUTED BENCHMARK

A solution can look smooth and still contain numerical error. Refine this mesh to see how a finite element temperature field approaches a known analytical solution.

GOVERNING EQUATION · 2D TRANSIENT HEAT CONDUCTION

\frac{\partial T}{\partial t} = α (\frac{\partial^{2} T}{\partial x^{2}} + \frac{\partial^{2} T}{\partial y^{2}})

T is the temperature rise above a uniform reference temperature, in kelvin (K).

A simple model with an exact answer

A homogeneous, isotropic rectangle cools by diffusion. All four edges are held at zero temperature rise. There are no internal heat sources.

We compare every mesh at the same time, t = 1 s, using the same color scale.

Domain: 0 ≤ x ≤ 1 m; 0 ≤ y ≤ 0.75 m
Diffusivity: α = 0.01 m²/s, constant
Initial condition: T(x,y,0) = 100 sin(πx/Lₓ) sin(πy/Lᵧ) K
Boundary condition: T = 0 K on all four edges, for t ≥ 0
Elements: Continuous bilinear quadrilaterals (Q1)

Analytical solution & numerical method

The separable sine-mode solution is:

T_{exact} = 100 \sin (\frac{π x}{L_{x}}) \sin (\frac{π y}{L_{y}}) e^{- α π^{2} (L_{x}^{- 2} + L_{y}^{- 2}) t} K.

Offline, scikit-fem assembles the consistent mass and stiffness matrices. The system MṪ + αKT = 0 uses nodal interpolation of the initial condition and Crank–Nicolson integration with Δt = 0.0025 s (400 steps). Boundary degrees of freedom are eliminated.

Relative L² error is integrated over the full Q1 field using quadrature order 8, checked against order 10. Halving Δt changes the field by less than 8.4 × 10⁻⁹ in relative L² norm—small compared with the plotted errors. The metric includes the remaining time-integration error; spatial discretization dominates here.

Only saved nodal values and metrics are loaded. The browser interpolates these values for display; it does not solve the heat equation.

Loading the local benchmark data…

What does the error measure?

The relative L² error measures the difference from the exact temperature field over the entire rectangle at t = 1 s. It is dimensionless; the readout shows it as a percentage.

E = \frac{{‖ T_{h} - T_{exact} ‖}_{L^{2}}}{{‖ T_{exact} ‖}_{L^{2}}} = \sqrt{\frac{\int_{Ω} {(T_{h} - T_{exact})}^{2} d A}{\int_{Ω} T_{exact}^{2} d A}}

Tₕ is the finite element field; Ω is the rectangle. The norm is evaluated by area integration, not just by comparing nodal values. Raw E = 0.01 means 1% relative error.

Refinement is a check, not a guarantee.

Discretization error & convergence

Finite elements approximate a continuous field. On this smooth benchmark, halving the element widths reduces error by about a factor of four, consistent with second-order L² convergence. Other models can converge differently.

Residuals answer different questions

An algebraic solver residual is the defect in the assembled discrete equations. A strong-form PDE residual is the defect when a numerical field is substituted into the differential equation. Neither is the error plotted here. Residual norms need their own definition and scaling and need not decrease monotonically across meshes.

Verification ≠ validation

Verification checks whether the equations are solved correctly; mesh convergence is one part. Validation compares model predictions with physical evidence. A fine mesh cannot fix incorrect physics, material data, or boundary conditions.

Download benchmark data Offline generation script