1. Errors and residuals

A solver can finish before your question is answered

Imagine predicting the temperature at the middle of a rod. The linear solver reports a tiny residual. That tells you something about a system of algebraic equations; it does not yet say whether the mesh resolves the temperature or whether the conductivity describes the material.

We will distinguish these questions using one matrix and one executable Poisson experiment. You need matrix multiplication and the idea of a vector length. The heat-transfer book supplies the physical interpretation if you want one.

Follow the approximations

Let $u$ solve the continuous mathematical problem, $u_h$ solve its discretized equations exactly, and $\widetilde u_h$ denote the field reconstructed from the computed coefficients. Then

$$u-\widetilde u_h=(u-u_h)+(u_h-\widetilde u_h).$$

The first term is spatial approximation error. The second includes incomplete algebraic solution and rounding. Neither term measures the difference between the mathematical model and the physical world.

Source of discrepancy	A controlled investigation
Physical assumptions	Compare a neglected effect with the retained balance
Uncertain input	Vary that input within a justified range
Spatial approximation	Refine the mesh with the model and solver policy fixed
Temporal approximation	Change time accuracy while keeping the spatial problem fixed
Algebraic iteration	Tighten the stopping tolerance on the same mesh
Floating-point arithmetic	Examine scaling and sensitivity; refine with care

These are different experiments. Refining a mesh cannot fix the wrong heat source, and tightening a solver cannot add missing spatial degrees of freedom.

What a residual measures

For $A u_h=b$, define the residual of a computed coefficient vector $v$ by

$$r=b-Av,\qquad e=u_h-v.$$

Subtraction gives $Ae=r$, hence $e=A^{-1}r$ when $A$ is invertible. Taking compatible norms yields

$$\frac{\lVert e\rVert}{\lVert u_h\rVert} \leq \kappa(A)\frac{\lVert r\rVert}{\lVert b\rVert}, \qquad \kappa(A)=\lVert A\rVert\lVert A^{-1}\rVert,$$

for nonzero $b$. The condition number describes how strongly residuals can be amplified into coefficient errors. Scaling and stopping criteria therefore matter together; see Barrett and coauthors’ discussion of stopping criteria.

A small residual with a large error

Take

$$A=\begin{pmatrix}1&0\\0&10^{-8}\end{pmatrix}, \quad b=\begin{pmatrix}1\\10^{-8}\end{pmatrix}, \quad u_h=\begin{pmatrix}1\\1\end{pmatrix}, \quad v=\begin{pmatrix}1\\0\end{pmatrix}.$$

Then $r=(0,10^{-8})^T$, although one component of the solution is completely wrong. In the infinity norm, the relative residual is $10^{-8}$ and the relative error is $1$. Here $\kappa_\infty(A)=10^8$, so the bound is attained.

If the two unknowns represent quantities with different units, a raw vector norm may itself be meaningless. Choose physical reference scales first: for example, measure temperature error relative to a chosen temperature difference.

Compute a spatial comparison

From the eqiora-source checkout created in Get started, run:

bash

cargo run -p eqiora-numerics --example poisson_convergence

The comparison program solves $-u''=\pi^2\sin(\pi x)$ with zero endpoint values on the unit interval. Before inspecting its table, differentiate $u=\sin(\pi x)$ twice and check the boundary values. That gives you the answer without relying on the program.

Read the fem_l2 and fvm_l2 columns as field errors, and the two balance columns as conservation discrepancies. Their magnitudes need not be similar: a method can balance its source and boundary flux almost exactly while its field remains visibly approximate. The following lessons derive both methods and explain precisely which field each error measures.

The program assembles and solves the system itself. Open its source to see where the model is loaded and where the mesh counts enter the comparison.

Estimate spatial error from an accepted result

The development checkout provides result.discretization_error_indicators(field, quadrature_points=3) for a steady scalar P1 triangular diffusion problem with constant positive diffusion and homogeneous essential boundaries. It returns an immutable array in the Result mesh’s cell order and its physical Dimension. The released 0.1.3 wheel predates this method.

For $-a\Delta u=f$ with constant $a>0$, each triangle receives

$$\eta_K^2=h_K^2\int_K (f/a)^2\,dx +\frac12\sum_{E\subset\partial K,\ E\text{ interior}} h_E\int_E [\nabla u_h\cdot n]^2\,ds.$$

Here $h_K$ is the longest triangle edge and $h_E$ is the edge length. Each interior-edge contribution is shared between its two triangles. The global estimate is $(\sum_K\eta_K^2)^{1/2}$ and estimates the H1 error seminorm; in two spatial dimensions its units are the Field’s units. The source-squared integral uses the explicitly selected Duffy quadrature, so the caller must establish its integration accuracy. See the residual estimator derivation for the element residual, normal jump and shape-regularity assumptions.

These are estimates, not certified upper bounds with a known constant. They do not include model discrepancy, quadrature error or an independently bounded algebraic error. Other boundary laws, variable diffusion, reaction, transport, transient problems and other Spaces reject rather than silently omit their residual terms. Mesh marking, local refinement and field transfer are not yet provided by this query.

The mesh error estimation notebook defines the Geometry, equations, mesh and solver in its own cells. It compares three uniform refinements with exact manufactured-solution errors, then saves and reopens the Plans and Results. Use the development checkout; hosted Colab execution has not been verified.

A misleading improvement

Suppose one run uses 16 cells and a relative solver tolerance of $10^{-3}$; another uses 32 cells and $10^{-10}$. A smaller field error in the second run does not isolate spatial convergence: two controls changed. First compare tolerances on the 16-cell problem, then compare meshes using an adequately tight common tolerance.

Conversely, a flat error curve under mesh refinement may mean algebraic error has become dominant. It can also mean a wrong comparison function or norm. Identify which hypothesis a follow-up run would distinguish.

Exercises

1. Replace $10^{-8}$ by $10^{-4}$ in the matrix example. Recompute the residual, condition number, and error without a solver.
2. Let $A=I$. Explain why residual and coefficient error now agree, and why this still says nothing about spatial discretization error.
3. In the Poisson output, find a row with a small balance discrepancy and a larger field error. Explain how both can be correct.
4. Design three runs that distinguish a solver-tolerance plateau from spatial error. State what stays fixed in each comparison.

Book map · Next: Time integration