# Which way is the right way to compute the integrals in finite element methods?

Finite element methods involve integrals of functions that are not polynomials, and these integrals must be computed numerically.

For example, suppose that $$f$$ is the right-hand side of a Poisson problem and we use Lagrange elements of degree $$k$$. The right-hand side of the discrete linear system is computed from integrals

$$\int_T f(x) p(x)$$

where $$p(x)$$ is a polynomial of degree at most $$k$$ over the element $$T$$. What is the right way of computing that integral?

1. Use numerical quadrature for the integral. Considering the case $$f = constant$$ suggests that the quadrature formula should exactness at least $$k$$.

2. Replace the function $$f(x)p(x)$$ by a polynomial interpolant $$\Pi[fp]$$ and compute $$\int_T \Pi[fp](x)$$ exactly. Again, the interpolant should have at least degree $$k$$.

3. Replace the function $$f(x)$$ by a polynomial interpolant $$\Pi[f](x)$$ and compute $$\int_T \Pi[f](x) p(x)$$ exactly.

Both quadrature rules and interpolation require point evaluations of $$f$$. However, finite element textbooks always review quadrature rules. What are the considerations here?

• You must also consider that curved elements introduce a Jacobian determinant when performing the transformation from the reference element (used to define integration rules). For linear elements this is constant and relates to the element volume. Otherwise, the Jacobian contributes to the polynomial order of the integrand. Jul 25, 2022 at 13:30

Quadrature and replacing the integrand by a polynomial are identical. That is how quadrature rules are derived.

To see why this is true remember that quadrature evaluates the integrand at only a select number of points. When you do quadrature, you cannot distinguish between function $$f(x)p(x)$$ and $$\Pi[fp](x)$$ (where $$\Pi[fp]$$ is the function that interpolates $$fp$$ and $$\Pi[fp](x)$$ is that function evaluated at $$x$$) in the sense that $$\int_a^b f(x) p(x) \, dx = \int_a^b \Pi[fp](x) \, dx$$ if you approximate both sides by quadrature. That's because during quadrature, you evaluate the integrands only at the interpolation points, and at these points the two integrands are the same by definition.

As a consequence, your approaches 1 and 2 are actually the same. Approach 3 is the same as the first two if you choose your quadrature points and weights appropriately.

• Thanks. On an interval, Gauss or Gauss-Labatto seems a reasonable choice. But there doesn't seem to be a tight correspondence between interpolation and quadrature on triangles and tetrahedra? Maybe even quads? For example, what is the interpolation corresponding to collapsed Gauss-quadrature on triangle? Perhaps you could comment on that. Jul 27, 2022 at 12:04
• @shuhalo The construction of quadrature points and weights works the same way on triangles and quadrilaterals. You interpolate by a polynomial and then you integrate the polynomial exactly. Jul 27, 2022 at 22:02

Personally I typically use 1, and something similar to 2/3 though it'd be more correct to say that I decompose $$f$$ into a polynomial series of degree $$k$$ and then integrate that multiplied by $$p$$ exactly. I think there might be some way to prove some sort of equivalence with the interpolation methods you described.

I tend to find that 1 gives me the most flexibility and for a nodal FEM scheme where my nodes are located at the quadrature points your method will still converge at the expected rates. If I really need to I can interpolate to a higher order quadrature scheme's nodes and use the same fundamental implementation with higher order quadrature. This is how I tend to measure L-2 convergence since in that case you need a quadrature method of at least $$O(2k+1)$$, and I typically go a little bit higher just in case.

The other method I use of decomposing $$f$$ into a basis of degree $$k$$ tends to work well when doing something like an explicit time integration scheme because then you can pre-compute your inverse mass matrix and multiply it with the other matrix of coefficients to get a simple matrix-vector product to solve for your fields: $$\partial_t q + \partial_x f = 0\\ \int_T p \partial_t q + \int_{\partial T} (f n)^* p - \int_T f \partial_x p = 0\\ \partial_t q_h = M^{-1} D f_h - M^{-1} L (f n)_h^*\\ M_{ij} = \int_T p_i(x) v_j(x)\\ D_{ij} = \int_T p_i(x) \partial_x v_j(x)\\ L_{ij} = \int_{\partial T} p_i(x) v_j(x)$$ You can of course do something similar with a second-order spatial derivative.

If I want to measure L-2 convergence I still need to perform higher-order quadrature, though.

• Thanks. Is there a textbook or an article that explains that? Jul 25, 2022 at 8:28
• The text I'm familiar with is Nodal Discontinuous Galerkin Methods by Hesthaven, he mainly discusses using the second method since explicit time stepping is used a lot in hyperbolic problems. He also presents the mathematical relationship between the various methods. Jul 25, 2022 at 16:28