Excursus - FEM - Gauss Quadrature + Example


#1

Numerical Integration

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 apriori.
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.

\underline{\text{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 (9) 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 you have to keep in mind is that

$$ \bbox[5px,border:2px solid red]
{
\text{for n int. points, Gauss integration provides 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 ---------- As an example see the following document: [Example 1](http://www.math.usm.edu/lambers/mat460/fall09/lecture31.pdf) If you want me to give you homework for this specific topic please let me know :) ---------- **Sources:** • Introduction to the Finite Element Method - Niels Ottosen & Hans Petersson • A First Course in Finite Elements by Jacob Fish and Ted Belytschko • Lectures of KIT Karlsruhe (Script: Einführung in die Finite Elemente Methode Sommersemester 2014) If you find any mistakes or have wishes for the next chapters please let me know.

#2

Example

I = \int\limits_{2}^5 (x^3 + x^2)\;dx \qquad 2m - 1 = 3 \Rightarrow m = 2

m = 2 \rightarrow two-point integration \rightarrow integral can be evaluated exactly

Formula (10) for (W_1,\xi_1) and (W_2,\xi_2):

\begin{bmatrix} 1 & 1 \\ \xi_1 & \xi_2 \\ \xi_{1}^2 & \xi_{2}^2 \\ \xi_{1}^3 & \xi_{2}^3 \end{bmatrix} \begin{bmatrix} W_1 \\ W_2 \\ \end{bmatrix} = \begin{bmatrix} 2 \\ 0 \\ \frac{2}{3} \\ 0 \\ \end{bmatrix} \Rightarrow \begin{gather} W_1 = W_2 = 1\\ \xi_1 = -\frac{1}{\sqrt{3}}\xi_2 = \frac{1}{\sqrt{3}} \end{gather}

By symmetry W_1 = W_2 and \xi_1 = -\xi_2. The first equation gives us the weights

1 \cdot W_1 + 1 \cdot W_2 = 2 \Rightarrow W_1 = W_2 = 1

as mentioned above. The third equation gives us the integration points.

We now use formula (2) with a = 2 and b = 5 to express x in terms of \xi:

x = \frac{1}{2}(a+b) + \frac{1}{2}\xi(b - a) = 3.5 + 1.5\xi \\ f(\xi) = (3.5 + 1.5\xi)^3 + (3.5 + 1.5\xi)^2

Using formula (5), the integral becomes:

$$\begin{align}
I& = J\hat{I} = \frac{l}{2} \int\limits_{-1}^1 ((3.5 + 1.5\xi)^3 + (3.5 + 1.5\xi)^2);d\xi \
& = \frac{l}{2}W_1((3.5 + 1.5\xi_1)^3 + (3.5 + 1.5\xi_1)^2) + \frac{l}{2}W_2((3.5 + 1.5\xi_2)^3 + (3.5 + 1.5\xi_2)^2) \
& = 37.818 + 151.432 = 191.25.
\end{align}

The result of the Gauss quadrature is exactly the same as if you would calculate the integral analytically.

\int\limits_{2}^5 (x^3 + x^2);dx = \left(\frac{x^4}{4} + \frac{x^3}{3}\right)_2^5 = 197.917 - 6.667 = 191.25

However in order not to calculate every $W_i$ and $\xi_i$ I will show you a table that gives you the weights and location depending on the number of Gauss points you chose. <img src="/forum/uploads/default/original/2X/9/94f059a5293592133f29670d749f2d7579cea427.gif" width="546" height="341"> ---------- **Source:** • A First Course in Finite Elements by Jacob Fish and Ted Belytschko • [Gauss Quadrature Table](http://www.unioviedo.es/compnum/laboratorios_web/Laborat10_inte/puntos_pesos_gauss.gif)

#3

@jousefm, perfect for linking to your blog post :smile: Great work!


#4

looks great!