Numerical methods. Note 3
[~Home~] Integration II. Numerical Methods. Note~ « 3 »

Adaptive algorithms are built on pairs of quadrature rules, a higher order rule (fx open-4) and a lower order rule (fx open-2). The higher order rule is used to compute the best approximation to an integral over a given interval. The difference between the results of the higher order rule and the lower order rule gives an estimate of the error in the approximation. If the error is larger than the required accuracy the interval is sub-divided into (two) smaller sub-intervals and the procedure applies recursively to the sub-intervals. The reuse of the function evaluations made at the previous step is very important for the efficiency of the algorithm. The equally-spaced abscissas naturally provide such a reuse.

Gaussian quadratures seek to obtain the best approximation to the integral ab w(x)f(x)dx  ≈  ∑i=1..n wi f(xi) with optimally chosen abscissas xi and weights wi. The optimal abscissas happen to be the roots of polynomials which are orthogonal on (a,b) with the weight-function w(x). For w(x) ≡ 1 these are the Legendre polynomials. Because the number of free parameters is 2n (n abscissas and n weights) the Gaussian quadratures are accurate up to 2n-1 order compared to only n-1 order for equally spaced abscissas. Unfortunately the optimal points can not be reused at the next iteration in an adaptive algorithm.

Gauss-Kronrod quadratures represent a compromise between equally spaced abscissas and optimal abscissas. They reuse n points from the previous iteration (n weights as free parameters) and add m optimal points (m abscissas and m weights as free parameters). Thus the accuracy of the method is n+2m-1. There are several special variants of these quadratures fit for particular types of the integrands.

QUADPACK is (probably) the best integration library, available at netlib.org. It is also compiled into GNU scientific library which is installed on our server, gunilo. An example of how to use it is here.

Problems

  1. Implement an adaptive integration routine with open 2-point and open 4-point formulas along the following lines:
    function [J,err] = adapt(fun,a,b,acc,eps)
    	[J,err] = open24(fun,a,b) ; tolerance = acc+eps*abs(J)
    	if ( err > tolerance ) then
    		[J1,err1] = adapt(fun,a,(a+b)/2,acc/2,eps)
    		[J2,err2] = adapt(fun,(a+b)/2,b,acc/2,eps)
    		J=J1+J2 ; err = err1 + err2
    
    One can use recursive functions in C, Fortran90 (here is an example) and Octave, but not in Fortran77. Reuse points !
  2. Calculate -2020 dx/1+x2 with acc=eps=0.0001 and compare the number of points (integrand evaluations) with DUMBINT.
  3. Non-obligatory. Compare your integrator with one of the QAG* routines from the QUADPACK library. The fortran version, qag.f, is here while the c version "gsl_integration_qag" is included in GSL (here is an example). Octave/Matlab can compare with a build-in integration routine.
    For example, calculate 01  ln(x)/√(x) dx = -4.

"Copyleft" © 2004 D.V.Fedorov (fedorov at phys dot au dot dk)