Numerical methods. Note 2
[ Home ] Numerical Integration Numerical Methods. Note  « 2 »
NB: Ingen undervisning i paaske uge!

Numerical integration rules (quadratures) seek to obtain an approximation

∫  0h f(x)dx  ≈  ∑i=1..n wi f(xi).

Newton-Cotes quadratures use equally spaced abscissas, e.g. xi=a+i/n+1h. The n free parameters w are then chosen such that the quadrature integrates exactly n basis functions, e.g. [1,x,x2,...,xn-1]. Examples of open rules:

2-point: x=[1/3h , 2/3h] , w=[h/2 , h/2];
4-point: x=[1/6h, 1/3h, 2/3h, 5/6h], w=[h/3, h/6, h/6, h/3].

Gaussian quadratures use optimally chosen weights w and also abscissas x. The optimal abscissas happen to be the roots of polynomials which are orthogonal on (a,b) with the weight-function w(x). Since 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.

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 approximation to the integral. The difference between 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.

Gauss-Kronrod quadratures represent a compromise between equally spaced abscissas and optimal abscissas: n points are reused from the previous iteration (n weights as free parameters) and then m optimal points are added (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, genryc. An example of how to use it is here.

Problems

  1. Implement an open integrator open24 wich calculates 2- and 4-point quadratures and returns the estimates of the integral and the error:
    function open24(f,a,b){ h=b-a; f1=f(a+h*1/6); f2=f(a+h*2/6); f3=f(a+h*4/6); f4=f(a+h*5/6); i2=(f2+f3)/2*h; i4=(2*f1+f2+f3+2*f4)/6*h; i=i4; err=abs(i4-i2); return {i:i,err:err};}
  2. Implement an adaptive integrator which estimates the integral with a given absolute (acc) and relative (eps) accuracy:
    function adapt(f,a,b,acc,eps){ y=open24(f,a,b); tol = acc+eps*abs(y.i); if (y.err < tol) then return y; else {y1=adapt(f,a,(a+b)/2,acc/√[2],eps); y2=adapt(f,(a+b)/2,b,acc/√[2],eps); return {i:y1.i+y2.i, err:√[y1.err2 + y2.err2]}}
  3. Test your implementation on some simple integrals. Reuse points.
  4. Calculate ∫01   ln(x)/√[x] dx = -4 with acc=eps=0.001 and estimate the number of need integrand evaluations.
  5. 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).

Questions