previous up next SYLLABUS  Previous: 3.2 An engineer's formulation  Up: 3 FINITE ELEMENT METHOD  Next: 3.4 Implementation and solution


3.3 Numerical quadrature

Slide : [ Formula - Table - Tunable || VIDEO login]

Unless the coefficients of the differential equation (here $ u$ and $ D$ ) are singular and require a special treatment somewhere in the solution domain, the precision of the quadrature depends mainly on the choice of the basis and test functions: it is indeed sufficient to simply preserve the convergence rate that is guaranteed by the FEM discretization. Writing the numerical quadrature as a truncated sum

$\displaystyle \int_a^b f(y) dy$ $\displaystyle =$ $\displaystyle \frac{b-a}{2}\sum_{i=1}^n w_i f(y_i) + R_n \hspace{10mm}$ (1)
$\displaystyle y_i$ $\displaystyle =$ $\displaystyle \left(\frac{b-a}{2}\right) x_i + \left(\frac{b+a}{2}\right)$ (2)

table 3.3.0#tab.1 shows how a judicious choice of the abscissas $ x_i$ and weights $ w_i$ yields an approximation of $ \mathcal{O}(f^{(p)})$ that integrates exactly polynomials with a degree up to $ (p-1)$ .

Table 3.3.0#tab.1: Quadrature abscissas and weights for the interval $ x_i\in [-1;1]$
quadrature $ n$ terms $ \pm x_i$ $ w_i$ $ R_n\sim\mathcal{O}(f^{(p)})$
mid-point 2 0 1 p=2
trapezoidal 2 1 1 p=2
Gaussian 2 $ \sqrt{1/3}$ 1 p=4
Gaussian 3 0 8/9 p=6
$ \sqrt{3/5}$ 5/9
Gaussian 4 $ \sqrt{3/7 +\sqrt{120}/35}$ $ 1/2 -5/(3\sqrt{120})$ p=8
$ \sqrt{3/7 -\sqrt{120}/35}$ $ 1/2 +5/(3\sqrt{120})$


For example, a superposition of $ n=4$ terms (where the integrand is evaluated on the specified location $ y_i(x_i)\in[a;b]$ and weighted by the factor $ w_i$ ) is required using a Gaussian quadrature, when four evaluations of the integrand in every mesh interval are used to compute exact overlap integrals involving the product of (two) cubic test and basis functions. In the case of linear roof-top functions, the mid-point and trapezoidal rules are both equally precise and require only one evaluation per interval; although this is slightly more expensive with one extra evaluation per interval, both rules can nicely be combined into a so-called tunable integration [8]

$\displaystyle \int_a^b f(y) dy = (b-a)\left[\frac{p}{2}\left[f(a)+f(b)\right] + (1-p) f\left(\frac{a+b}{2}\right)\right] + R_2, \hspace{5mm}p\in[0;1]$ (3)

As will become clear below, this allows to gradually modify the piecewise linear FEM discretization (obtained with $ p=1/3$ ) into either an equivalent FD discretization (with $ p=1$ ) or a piecewise constant FEM approximation (with $ p=0$ ). Apart from the academic interest, this feature can be used to change the slope of the numerical convergence and even to stabilize an approximation that is marginally unstable.

Armed with a variety of quadrature rules, it is now possible to implement the linear FEM scheme sketched in sect.3.2 using a tunable integration (3.3#eq.3) to evaluate the matrix elements in (3.2#eq.5). Since roof-top functions reach only as far as their nearest neighbors, all the matrix elements $ \mathbf{A_{ij}}$ with $ \vert i-j\vert>1$ vanish, except for the overlap at the boundaries. Choosing a sequential numbering of the unknowns $ \{f_j^{t+\Delta t}\}, j=1,N$ , this results in a tri-diagonal matrix, with two extra elements in the upper-right and lower left corner taking care of the periodicity. It is a good exercise for the reader to verify that all the matrix elements can be written in terms of the integrals

$\displaystyle \int_{x_{j-1}}^{x_{j}}e_{j}e_{j}dx =\int_{x_{j}}^{x_{j+1}} e_{j}e_{j} dx = (1+p)\frac{\Delta x}{4},$    $\displaystyle \int_{x_{j}}^{x_{j+1}}\!\!e_{j}e_{j+1}dx =\int_{x_{j}}^{x_{j+1}}\!\!e_{j+1}e_{j} dx = (1-p)\frac{\Delta x}{4}$  
$\displaystyle \int_{x_{j}}^{x_{j+1}}e_{j}e_{j}^\prime dx = -\int_{x_{j}}^{x_{j+1}}e_{j}e_{j+1}^\prime dx = -\frac{1}{2},$    $\displaystyle \int_{x_{j}}^{x_{j+1}}e_{j}^\prime e_{j}^\prime dx = -\int_{x_{j}}^{x_{j+1}}e_{j}^\prime e_{j+1}^\prime dx = \frac{1}{\Delta x}$  

provided that the advection and diffusion coefficients are constant across the intervals as will be assumed here for pedagogical reasons. In a real code constructed for an inhomogeneous medium, each formula is replaced by $ n$ function calls returning the value of the integrand (where $ u(x_k)$ , $ D(x_k)$ are evaluated on quadrature points) and the weighted sum (3.3#eq.1) is simply accumulated in the corresponding matrix element.

SYLLABUS  Previous: 3.2 An engineer's formulation  Up: 3 FINITE ELEMENT METHOD  Next: 3.4 Implementation and solution

      
back up next contents bibliography Copyright © Lifelong-learners at 01:45:59, February 11th, 2012