Solving a linear PDE like $$ \Delta u = f \quad\text{on } \Omega,\\ n\cdot \nabla u = 0 \quad\text{on } \Gamma, $$ with Finite Elements usually goes like this:

  1. Create the discretization $Au=b$ via $$ \int \nabla u \cdot \nabla v_i = \int f v_i $$ (With the $u$ and $v_i$ from some finite-dimensional trial and test spaces; let's say they the spaces are equal.)

  2. Forget everything about the continuous problem and simply solve $Au=b$ with a method of your choosing, e.g., some Krylov method.

When solving the discretized equation system with an iterative solver, often some norm of the discrete residual is minimized; in the case of MINRES and GMRES, for example, the discrete 2-norm, $$ \|r|\|_2^2 = \|b - Au\|_2^2 = \sum_i \left(\int f v_i - \nabla u \cdot \nabla v_i\right)^2. $$


Does this norm have any correspondence in the non-discretized world, i.e., is there an entity in terms of $u$ and $f$ to which $\|r\|_2$ is an approximation?

  • 5
    $\begingroup$ Yes, the residual of the weak form of the PDE (i.e., the thing that acts on test functions), but it is a poor approximation since a) it doesn't scale correctly under mesh refinement (for which you could use a weighted 2-norm) and b) the 2-norm is not the appropriate norm in which to measure the weak residual. This is closely tied to the question of preconditioners; you might be interested especially in etna.mcs.kent.edu/volumes/2011-2020/vol41/… and arxiv.org/abs/1603.04475. $\endgroup$ Commented Apr 23, 2018 at 18:30

2 Answers 2



If you're using scalar products in FEM/FVM discretizations, use the mass-matrix scalar product, not $\ell_2.$


If you're solving FEM/FVM systems with Krylov methods, precondition with the inverse of the mass-matrix.

Some details

It's easy to think that the 2-norm (or in fact any scalar product) $\|x\|_2$ of a finite-something discretization is an approximation of the $L_2$-norm of the corresponding function. That is not the case.

Let $v_i$ be the basis functions of the finite-dimensional space $V_h$, and let the coefficients be denoted by Greek letters such that all functions in $V_h$ can be written as $v = \sum_i \alpha_i v_i$. Then the natural $L_2$ scalar product is $$ \langle u, v\rangle_2 = \int u v = \sum_i\sum_j \int (\alpha_i v_i) (\beta_j v_j) = \alpha^T M \beta = \langle\alpha, \beta\rangle_M $$ with $M$ being the mass matrix, $M_{ij} = \int v_i v_j$. This is the scalar product that Krylov methods of discretizations should be formulated in, not $\alpha^T \beta$ as it's usually done. This has multiple consequences:

  • Symmetry is no longer the determining propery of $A$, it's self-adjointness in $\langle\cdot,\cdot\rangle_M$: Instead of $A$, one needs to look at the system matrix $M^{-1}A$.

  • The residual is $r=M^{-1} F_h(u)$; calling $\tilde{r}=F_h(u)$ our "old" residual, one has $$ \langle r, r\rangle_M = \langle M^{-1} \tilde{r}, M^{-1} \tilde{r}\rangle_M = \tilde{r}^T M^{-1} \tilde{r}. $$ In the case of finite volumes, it's easy to see that this quantity actually approximates the continuous resisdual $\int F(u)^2$: The mass matrix is $M_ij = \delta_{i,j} |\Omega_i|$ (a diagonal matrix with the finite volumes). Since for small $\Omega_i$, $F(u)$ changes very little across $\Omega_i$ (assuming it's continuous), it can be taken out of the integral $$\begin{split} \tilde{r}^T M^{-1} \tilde{r} &= \sum_i |\Omega_i|^{-1} \left(\int F(u) v_i\right)^2\\ &= \sum_i |\Omega_i|^{-1} \left(\int_{\Omega_i} F(u)\right)^2\\ &\approx \sum_i |\Omega_i|^{-1} F(u)(x_i)^2 \left(\int_{\Omega_i} 1 \right)^2\\ &= \sum_i |\Omega_i| F(u)(x_i)^2\\ &\approx \int F(u)^2. \end{split}$$ For finite elements, it's harder to see. (Perhaps someone else can fill that in.)

The article http://etna.mcs.kent.edu/volumes/2011-2020/vol41/abstract.php?vol=41&pages=13-20 explains how replacing the scalar product is equivalent to preconditioning.


I am digressing slightly to address how the equation is discretized before going to the solution.

General methods for numerically solving PDEs:

  1. Use pseudospectral methods: Approximate the solution as a series expansion (often Fourier or Chebyshev). Most of these series form an orthonormal basis set. Discretize the PDE as a set of matrix equations. This is possible because the differentiation of series expansion is rather straightforward. Solve the matrix equation.

Here is an example using the strong form of the differential equation. $T_n(x)$ represents the Chebyshev polynomial. You could replace that with any other basis set of your choice.

$$ \nabla^2u = f \\ Let\ u = \sum_{i=1}^Na_iT_{i-1}(x) $$ Here, $a_i$ are the coefficients we need to find. Substituting the value of $u$ in the differential equation, we get

$$ \nabla^2\sum_{i=1}^Na_iT_{i-1}(x) = f \\ \sum_{i=1}^Na_i\nabla^2T_{i-1}(x) = f $$

Then, we discretize $x$ and $f$ on a set of collocation points ($0$ to $N$). This leads us to the matrix equation $$ L.a = f $$ where, $L$ is the discretization operator, which contains the second derivatives of the Chebyshev series. In essence, we are choosing a set of $N$ points, where we make the residual zero.

  1. Use Galerkin methods: In addition to approximating the solution as a series expansion, we project our solution onto another vector space, by multiplying the equation with a test function and integrating over the domain. For Galerkin methods, the trial function used to approximate the solution is the same as the test function used for projection. The use of variational form instead of the strong form of PDE is somewhat historical, and somewhat out of convenience. The main advantage is that by integrating and using the Green's theorem, you can reduce one order of the differential equation. Therefore, you do not have to calculate higher order derivatives. The matrix equation you would get is the same as the equation above.

The only difference between finite element (Galerkin) and pseudospectral methods is that the test function in the pseudospectral method is a Dirac delta function. Since you already have the equations for the finite element approximation in the question, I am not writing them here.

So the general matrix equation we get from any method is: $$ Lu = f $$

Remember that $u$ here is not the solution, but the coefficients of the solution approximated as series expansion.

The most direct way to solve this equation using iterative methods is to solve the time-dependent equation: $$ u_t = -Lu + f $$

This is a generalized diffusion equation. Here the time is actually the number of iterations, and we use the analogy of a diffusion equation to iteratively solve the matrix equation. The analogy of iteration = diffusion is very useful because all iterative methods "diffuse" the error until the solution settles down to a steady state value.

Any stationary, one step iterative method will have the form of $$ u_{n+1} = Gu_n + k\\ G = I - R^{-1}L \\ k = R^{-1}f $$ where $I$ is the identity matrix, and $R$ is some easily invertible square matrix (a.k.a preconditioner). The values of $G, k$ would vary depending on the iterative method used. In principle, matrix functions can always be evaluated by using an expansion in terms of the eigenfunctions of the matrix. For the respective eigenvectors and eigenvalues $\phi_i, \lambda_i$ of the square matrix $G$, the stationary, one step method above can be written as: $$ u_{n+1} = \sum_{j=1}^{N}\lambda_ju^n_j\phi_j + k $$

We define the error at the nth step as: $$ e_n = u - u_n $$

Substituting this back to the matrix equation, we get: $$ e_{n+1} \equiv Ge_n $$

The magnitude of the largest eigenvalue of the iteration matrix $G$ is called the spectral radius of the matrix. Optimizing an iteration is equivalent to minimizing the spectral radius. You can easily see how the spectral radius or the matrix $G$ is related to the coefficients of expansion of our original differential operator. It is related to the set of basis functions, the specific linear combination of which is our solution. What you are minimizing is the coefficient of our approximate solution.

EDIT: Source: Chapter 15 http://depts.washington.edu/ph506/Boyd.pdf

  • 2
    $\begingroup$ Welcome to SciComp.SE! This is an extensive answer, but how precisely does it address the question about the continuous (function space) interpretation of the residual? $\endgroup$ Commented Apr 24, 2018 at 19:39
  • $\begingroup$ Thanks for welcoming! The interpretation of the residual is the spectral radius. It can be visualized as how much the linear operator is transforming/stretching the solution. I have tried to show how it is related to the function. I don't know of a better way to put it. $\endgroup$ Commented Apr 24, 2018 at 20:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.