Numerical Integration

One of the most widely used methods of numerical integration is Gauss-Legendre quadrature. It posses very attractive property of to be exact on polynomials of degree up to $2n-1$ , while using only $n$ integrand evaluations ( $n$ -point quadrature).
The algorithm consists in approximation of initial definite integral by the sum of weighted integrand values sampled at special points called abscissas:

$\begin{matrix} \displaystyle \int_a^b f(x)\,\mathrm{d}x\approx\frac{b-a}{2}\sum_{i=1}^{n}{w_i\,f\left(\frac{b-a}{2}\xi_i+\frac{b+a}{2}\right)}&\\ &\\ \displaystyle w_i = \frac{2}{\left(1-\xi_i^2\right)\,\left[P'_n(\xi_i)\right]^2}&(n=1,2,\dots) \end{matrix}$

where values $\xi_i$ are $n$ zeroes of the $n^{th}$ -degree Legendre polynomial $P_n(\xi)$ .
Accuracy of the numerical integration depends significantly on precision of zeroes $\xi_i$ . There is no easy-to-use analytical expression for $\xi_i$ and usually they are computed numerically by root-finding algorithms.
Surprisingly popular mathematical references/software provide low-precision $\xi_i,\,w_i$ (only 4-6 correct digits) for limited $n\le 5$ . Obviously such situation doesn’t reflect contemporary computational capabilities and accuracy demands of many applications. For instance, commonly used floating-point type double of IEEE 754 standard is capable to store values with 16 digits precision (machine epsilon is about 1e-16 for 64 bits double).
This page aims to provide software libraries for calculation of high-precision abscissas $\xi_i$ and weights $w_i$ for any desired $n$ .

Gauss-Legendre Quadrature for C/C++

This open-source library implements numerical integration based on Gauss-Legendre quadrature of any order. Pre-calculated high-precision abscissas and weights are used for specific orders n=2,…, 20, 32, 64, 96, 100, 128, 256, 512, 1024. Values for all other n are generated on the fly with 1e-10 precision by fast root-finding algorithm.

Also library includes routine for numerical integration over 2D rectangle using product of two 1D Gaussian quadratures.
If you are looking for numerical integration over the unit disk (2D sphere) you might be interested in this page Cubature formulas for the unit disk. It contains derivation details and source code in C/C++.

Gauss-Legendre Quadrature for MATLAB

There are several routines for Gaussian quadrature implemented using Matlab programming language. All of them suffer from low performance and even some of them provide low accuracy results.
To alleviate the situation I am developing quadrature toolbox which has Matlab interface (m-files) and delegates all the heavy lifting to high-performance MEX library written in C/C++ (see previous section). So far it includes only Gauss-Legendre rule with many others to come.
Quadrature toolbox consists of three files:

GaussLegendre.m
quadlab.mexw32
quadlab.mexw64

It was created and tested using Matlab R2010b, both on 32 and 64 bits versions. Usage example:

>>GaussLegendre(@cos,-pi/2,pi/2,1024)
  ans =
       2.000000000000000

GaussLegendre() accepts Matlab function (inline or from m-file), integration boundaries $[a,b]$ , order of quadrature (1024 in example) and desired precision for nodes as optional parameter. Also it returns vectors of abscissas and weights if such outputs are specified.
Performance comparison with the fastest Matlab-language routine (x-quadrature order, y-execution time, green-quadrature toolbox, red-competitor, lower is better):

Besides being fast library usually provides results with better accuracy thanks to pre-calculated high-precision abscissas and weights.

from http://www.holoborodko.com/pavel/numerical-methods/numerical-integration/

Gleaning on the Westward Journey

Search This Blog

Sunday, April 10, 2011

Gaussian Quadrature

Numerical Integration

Gauss-Legendre Quadrature for C/C++

Gauss-Legendre Quadrature for MATLAB

No comments:

Post a Comment