The correct way to compute errors

mathematics
1

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.

2

Yes. You can also define the quadrature rule using the metadata option in ufl.dx.

Again, then it leads to the expansion of the inner product.

You avoid round-off errors, at the cost of using the interpolated solution (and not the exact solution at quadrature points).

3

Ahhh so the main idea is to subtract nodal values, and the interpolation into the high-order space is just a necessary extra step we do to compensate the loss we get from the interpolation!
Thanks a lot, I appreciate your fast response.

4

How do I compute the error without interpolating the exact solution?

Say I have an exact solution defined like this:

# exact solution
class u_exact:
    def __init__(self):
        self.t = 0.0
        
    def eval(self, x):
        return np.sin(np.pi*x[0]) + 10*self.t

Can I avoid the interpolation in the high order polynomial space? Can I compute the error using the exact representation?

The reason I am asking is the Gibbs oscillation phenomenon.