Numerical Integration - Gauss Quadrature

In general, we cannot integrate the weak form of the FEM in closed form.


Or more precisely: The reason why we cannot integrate the polynomial trial functions analytically is that functions like E(x) or A(x) are not known beforehand.
To avoid this problem we need numerical integration. We approximate the real integral I with an approximation \tilde{I}.

In quadrature we analyse the integrand at fixed points \xi_i. The next step is to weight the calculated points wih \lambda_i and sum them up in order to get the approximation \tilde{I}.

\int\limits_{a}^b f(x)\;dx \approx \tilde{I} = \sum_{i=1}^m f(\xi_i) \lambda_i \tag{1}

The Gauss quadrature formula is always given in a parent domain [-1,1] to the physical domain [a,b]. Using the coordinate transformation

x = \frac{a + b}{2} + \xi\frac{b - a}{2} \tag{2}

and

\frac{dx}{d\xi} = \frac{b - a}{2} \Rightarrow dx = \frac{b - a}{2} d\xi \tag{3}

the integral can be mapped to the parent domain

I = \int\limits_{a}^b f(x)\;dx = \frac{b - a}{2} \int\limits_{-1}^1 f(\xi)\;d\xi \tag{4}

where the coefficient \frac{b - a}{2} is the so called Jacobian J.

I will now outline, how to setup the Gauss quadrature.
We approximate our integral by

\tilde{I} = \lambda_{1}f(\xi_{1}) + \lambda_{2}f(\xi_{2}) + \dotsc = \underbrace{\begin{bmatrix} \lambda_1 & \lambda_2 & \dotsc & \lambda_n \end{bmatrix}}_{\mathbf{\lambda}^\intercal} \underbrace{\begin{bmatrix} f(\xi_{1}) \\ f(\xi_{2}) \\ \vdots \\ f(\xi_{n}) \end{bmatrix}}_{\mathbb{f}} = \mathbb{\lambda}^\intercal \mathbb{f} \tag{5}

where \lambda_i are the weights as mentioned previously and \xi_i the points at which the integrand is evaluated.

The basic idea is simple:

Choose the weights and the integration points so that the possible polynomials are integrated exactly.

f(\xi) = \alpha_1 + \alpha_2\xi + \alpha_3\xi^2 + \dotsc = \underbrace{\begin{bmatrix} 1 & \xi & \xi2 & \dotsc \end{bmatrix}}_{\mathbf{p}} \underbrace{\begin{bmatrix} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \vdots \end{bmatrix}}_{\mathbb{\alpha}} = \mathbb{p(\xi)}\alpha \tag{6}

The next step is to express the values of the coefficients \alpha_i in terms of the function f(\xi) at the integration points:

\underbrace{\begin{bmatrix} f(\xi_1) \\ f(\xi_1) \\ \vdots \\ f(\xi_n) \end{bmatrix}}_{\mathbb{f}} = \underbrace{\begin{bmatrix} 1 & \xi_1 & \xi_{1}^2 & \cdots \\ 1 & \xi_2 & \xi_{2}^2 & \cdots \\ \vdots & \vdots & \vdots & \vdots \\ 1 & \xi_n & \xi_{n}^2 & \cdots \\ \end{bmatrix}}_{\mathbb{M}} \underbrace{\begin{bmatrix} \alpha_1 \\ \alpha_2 \\ \vdots \\ \alpha_n \end{bmatrix}}_{\mathbb{\alpha}} \tag{7}

Based on (5) and (7), the integral \hat{I} will be written as

\hat{I} = {\mathbf{W}}^\intercal \mathbf{M} \alpha \tag{8}

To determine how the weights and quadrature points look like, we integrate the polynomial f(\xi):

\tilde{I} = \int\limits_{-1}^1 f(\xi)\;d\xi = \int\limits_{-1}^1 \begin{bmatrix} 1 & \xi & \xi^2 & \xi^3 & \cdots \end{bmatrix} \begin{bmatrix} \alpha_1 \\ \alpha_2 \\ \vdots \\ \alpha_n \end{bmatrix} d\xi = {\begin{bmatrix} \xi & \frac{\xi^2}{2} & \frac{\xi^3}{3} & \frac{\xi^4}{4} & \cdots \end{bmatrix}}_{-1}^{1} \mathbf{\alpha} = \underbrace{\begin{bmatrix} 2 & 0 & \frac{2}{3} & 0 & \cdots \end{bmatrix}}_{\hat{\mathbf{P}}} \mathbf{\alpha} = \mathbf{\hat{P}} \mathbf{\alpha} \tag{9}

The weights and quadrature points are selected in a way that \hat{I} in the equation we just derived equals (8) so that we get the exact integral for a polynomial of a given order. This yields

\mathbf{W}^\intercal \mathbf{M} \alpha = \mathbf{\hat{P}} \mathbf{\alpha} \Rightarrow \boxed{\mathbf{M}^\intercal \mathbf{W} = \mathbf{\hat{P}}^\intercal} \tag{10}

The most important thing we have to keep in mind is that

$$ \bbox[5px,border:2px solid black]
{
\text{for n int. points, Gauss int. provides an exact integration of polynomials of the order 2n - 1.}
}

## Properties of Gauss Quadrature \+ weights always positive \+ optimal order of integration is $p = 2m -1$ m: = number of Gauss points \+ very good if you want a high accuracy with a low number of Gauss points (if we do not require an equidistant distribution) • Gauss points always distributed symmetrically, i.e. the mid of our interval is for an odd number of Gauss points always a point we can evaluate ---------- **For an example, please have a look at the following page:** **[Gauss Quadrature Example](https://www.simscale.com/forum/t/excursus-fem-gauss-quadrature-example/29535/2)** ---------- Also see our ​**[SimScale SimWiki](https://www.simscale.com/docs/content/simwiki.html)** for more about other interesting simulation related questions. ---------- ### **Literature References:** * Introduction to the Finite Element Method - Niels Ottosen & Hans Petersson * A First Course in Finite Elements by Jacob Fish and Ted Belytschko