.. !split

Numerical integration  (1)
==========================

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,

.. math::
        
        \int_{-1}^{1} g(X)dX \approx \sum_{j=0}^M w_j g(\bar X_j),
        

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

The very simplest method is the *Midpoint rule*,

.. math::
        
        \int_{-1}^{1} g(X)dX \approx 2g(0),\quad \bar X_0=0,\ w_0=2,
        

which integrates linear functions exactly.

.. _fem:approx:fe:numint1:

Newton-Cotes rules
------------------

.. index::
   single: numerical integration; Midpoint rule

.. index::
   single: numerical integration; Trapezoidal rule

.. index::
   single: numerical integration; Simpson's rule

.. index::
   single: numerical integration; Newton-Cotes formulas

.. index:: Midpoint rule

.. index:: Trapezoidal rule

.. index:: Simpson's rule

.. index:: Newton-Cotes rules

The `Newton-Cotes <http://en.wikipedia.org/wiki/Newton%E2%80%93Cotes_formulas>`__
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*,

.. _Eq:fem:approx:fe:numint1:trapez:

.. math::
   :label: fem:approx:fe:numint1:trapez
        
        \int_{-1}^{1} g(X)dX \approx g(-1) + g(1),\quad \bar X_0=-1,\ \bar X_1=1,\ w_0=w_1=1,
         
        

and *Simpson's rule*,

.. math::
        
        \int_{-1}^{1} g(X)dX \approx \frac{1}{3}\left(g(-1) + 4g(0)
        + g(1)\right),
        

where

.. math::
        
        \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 <http://tinyurl.com/nm5587k/fem/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 :math:`\int_{a}^{b}f(x)dx`.

Gauss-Legendre rules with optimized points
------------------------------------------

.. index:: Gauss-Legendre quadrature

More accurate rules, for a given :math:`M`, arise if the location of the
integration points are optimized for polynomial integrands.  The
`Gauss-Legendre rules <http://en.wikipedia.org/wiki/Gaussian_quadrature>`__ (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

.. math::
        
        M=1:\quad \bar X_0=-\frac{1}{\sqrt{3}},\ 
        \bar X_1=\frac{1}{\sqrt{3}},\ w_0=w_1=1
        

.. math::
          
        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 :math:`M`-point Gauss-Legendre rule integrates a polynomial
of degree :math:`2M+1` exactly.
The code ``numint.py`` contains a large collection of Gauss-Legendre rules.

.. 2DO

.. Newton-Cotes: bedre med det som overskrift over

.. Later:

.. lumped mass via num int; example or exercise

.. non-uniform meshes

.. hand-calculation, extend software?

.. example: half a Gaussian hat with one fine-grid area and a coarse-grid area

.. adaptivity