In this tutorial, we learn that if we are interested in finding the error of a numerical solution u_h to a given problem with exact solution u_{ex}, we **should not** compute directly the integral

\displaystyle \epsilon := \int_\Omega | u_h - u_{ex} |^2

which corresponds to

```
V = FunctionSpace(mesh, ("Lagrange", degree))
error = form((uh - u_ex)**2 * ufl.dx)
E = np.sqrt(comm.allreduce(assemble_scalar(error), MPI.SUM))
```

**Question 1**: I assume, that `form()`

is the function that calls the quadrature function in the background, which finds an approximation (depending on the quadrature rule and the degree of u_h and u_{ex}) to the integral above, right?

The reason is that after expanding the square root in the integral, the compute tries to integrate

\displaystyle u_{ex}^2 + u_h^2 - 2u_{ex}u_h

where the terms u_{ex}^2 + u_h^2 and 2u_{ex}u_h are of a similar size, if u_{ex} and u_h are of similar size (which is the case when the error is small) - leading to round-off errors.

**Instead of that we should** interpolate both u_h and u_{ex} into a higher order space of polynomials.

```
W = FunctionSpace(mesh, (family, degree + degree_raise))
u_W = Function(W)
u_W.interpolate(uh)
u_ex_W = Function(W)
u_ex_W.interpolate(u_ex)
```

and then we subtract the nodal values (the expansion coefficients of the FEM functions) of the two functions and compute the error as before:

```
e_W = Function(W)
e_W.x.array[:] = u_W.x.array - u_ex_W.x.array
error = form(ufl.inner(e_W, e_W) * ufl.dx)
E = np.sqrt(comm.allreduce(assemble_scalar(error), MPI.SUM))
```

**Question 2**: why do we subtract the arrays of nodal values instead of subtracting the `Function`

objects inside the `ufl.inner()`

like we did for the naive error computation?

**Question 3**: Why should this procedure give a more accurate error? Morally speaking, we are loosing information by interpolating the exact solution in a polynomial space, right? I tried to work it out, but I can not see why this approach avoids round-off errors.