Numerical integration

Finite element codes usually apply numerical approximations to integrals. Since the integrands in the coefficient matrix often are (lower-order) polynomials, integration rules that can integrate polynomials exactly are popular.

The numerical integration rules can be expressed in a common form,

\[\tag{107} \int_{-1}^{1} g(X){\, \mathrm{d}X} \approx \sum_{j=0}^M w_j g(\bar X_j),\]

where \(\bar X_j\) are integration points and \(w_j\) are integration weights, \(j=0,\ldots,M\). Different rules correspond to different choices of points and weights.

The very simplest method is the Midpoint rule,

\[\tag{108} \int_{-1}^{1} g(X){\, \mathrm{d}X} \approx 2g(0),\quad \bar X_0=0,\ w_0=2,\]

which integrates linear functions exactly.

Newton-Cotes rules

The Newton-Cotes rules are based on a fixed uniform distribution of the integration points. The first two formulas in this family are the well-known Trapezoidal rule,

\[\tag{109} \int_{-1}^{1} g(X){\, \mathrm{d}X} \approx g(-1) + g(1),\quad \bar X_0=-1,\ \bar X_1=1,\ w_0=w_1=1,\]

and Simpson’s rule,

\[\tag{110} \int_{-1}^{1} g(X){\, \mathrm{d}X} \approx \frac{1}{3}\left(g(-1) + 4g(0) + g(1)\right),\]

where

\[\tag{111} \bar X_0=-1,\ \bar X_1=0,\ \bar X_2=1,\ w_0=w_2=\frac{1}{3},\ w_1=\frac{4}{3}{\thinspace .}\]

Newton-Cotes rules up to five points is supported in the module file numint.py.

For higher accuracy one can divide the reference cell into a set of subintervals and use the rules above on each subinterval. This approach results in composite rules, well-known from basic introductions to numerical integration of \(\int_{a}^{b}f(x){\, \mathrm{d}x}\).

Gauss-Legendre rules with optimized points

More accurate rules, for a given \(M\), arise if the location of the integration points are optimized for polynomial integrands. The Gauss-Legendre rules (also known as Gauss-Legendre quadrature or Gaussian quadrature) constitute one such class of integration methods. Two widely applied Gauss-Legendre rules in this family have the choice

\[\tag{112} M=1:\quad \bar X_0=-\frac{1}{\sqrt{3}},\ \bar X_1=\frac{1}{\sqrt{3}},\ w_0=w_1=1\]

\[\tag{113} M=2:\quad \bar X_0=-\sqrt{\frac{3}{{5}}},\ \bar X_0=0,\ \bar X_2= \sqrt{\frac{3}{{5}}},\ w_0=w_2=\frac{5}{9},\ w_1=\frac{8}{9}{\thinspace .}\]

These rules integrate 3rd and 5th degree polynomials exactly. In general, an \(M\)-point Gauss-Legendre rule integrates a polynomial of degree \(2M+1\) exactly. The code numint.py contains a large collection of Gauss-Legendre rules.