$L^2$-error in FEM: how to compute integral over reference element?

Question

I have the following problem. The domain is $(0,1)$ and we consider a uniform triangulation on $\hat{\Omega}$ with elements $K_i = [i/N,(i+1)/N]$ and $X_h^1$ the linear finite element space. I wrote the following affine map to the reference element $\hat{K} = [0,1]$:

$$ x = x_i + \gamma (x_{i+1} - x_i),$$

the following basis function $\phi(x)$:

\begin{cases} \frac{x-x_{i-1}}{x_i-x_{i-1}} & x_{i-1} \leq x \leq x_i \\ \frac{x_{i+1}-x}{x_{i+1}-x_i} & x_i \leq x \leq x_{i+1} \\ 0, \end{cases}

end the following expression for the interpolant:

$$ \Pi_h^1 f(x) = \sum_{i=0}^{N+1} f(x_i) \phi_i(x). $$

I interpolated the function $f(x) = cos(2 \pi x)$ with a small program in Python and now I have a problem with the $L^2$-error defined as:

$$ ||f - \Pi_h^1 f(x)||^2 = \int_0^1 \left (f - \Pi_h^1 (f) \right )^2 dx$$

which can be decomposed into the sum of the integrals over each element of the mesh and each integral is computed on the reference element using the affine map. I have a problem with the part. I did something like this:

$$ \sum_{i=0}^{N_e} \left (\int_0^1 (f_e - \Pi_h^1 f_e) \right )^2 dx, $$

and even if this is correct, I do not know how to compute each integral on the reference element. Can somebody give me an hint, especially on how to compute $f$ on the reference element? In particular, I tried to wrote the integrand in the following way:

def func_ref(z, x, a, b):
    return np.power((np.cos(x[a]) - np.cos(x[a]) * (1 - z) +
    np.cos(x[b]) * z) * (x[b] - x[a]), 2)

where z is the variable.

Thank you.

Answer 1

You can write

$$\int_{0}^1 (f - \Pi_h^1(f))^2 dx = \sum_{i=0}^{N_e} \int_{e} (f - \Pi_h^1(f))^2 dx$$

where $e = [x_k, x_{k+1}]$ is an interval (triangle as you say).

Method A (simple, no reference element)

Now on each $e$, you know what is $\Pi_h^1(f)$:

$$\Pi_h^1(f)(x) = f(x_k)\frac{x_{k+1}-x}{x_{k+1}-x_{k}} + f(x_{k+1})\frac{x-x_k}{x_{k+1}-x_{k}}$$

and you can integrate using you favorite numerical quadrature.

Method B (reference element)

Using the reference element will lead you to this procedure:

• You start by defining the quadrature (nodes and weights) on the reference element.
• Now you want to evaluate $(f - \Pi_h^1(f))^2$ on the reference element. $\Pi_h^1(f)$ is quite easy, because it is a combination of basis functions (whose values you know in the reference element). For $f$, you have no choice but to map the quadrature point to the element $e$ (because this is where $f$ is defined).
• You gain one term which is the determinant of the inverse of the jacobian of the mapping (because of the coordinate change).

Finally, you can realize that you have evaluated the function in the same points, and the additional term from the jacobian of the mapping in method (B) is the one you use to adapt the weights of the quadrature rule in method (A). So, actually, it's just two ways of seeing exactly the same stuff.