Numerical methods. Note 6
Home ] Integration I Numerical Methods. Note  « 6 »

The classical n-point formula gives an estimate of an integral 0h f(x)dx ≈ In=∑i=1:nwif(xi) using the values of the integrand at n-points equally spaced within the interval (0,h). The weights wi are chosen such that the formula integrates precisely n polynomials x0,..,xn-1. Closed (Newton-Cotes) formulae use the end-points of the interval: Open formulae use only the inner points of the interval:

One can choose the weights wi (for open formulae) such that the formula integrates precisely diverging integrals, say, 1/√(x).

Due to round-off errors n cannot be too large, say n<10 for single and n<20 for double precision. Therefore if the accuracy is not good enough the interval is subdivided into smaller intervals and the formulae are applied to subintervals thus reducing the error.

The error can be estimated by comparing the results from formulae of different precision, e.g., open-2 and open-4.

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.

Problems

  1. Make a function OPEN24 that estimates the integral J=ab f(x)dx and the the integration error ERR using the 2-4 point open formulae:

    function [J,ERR] = OPEN24(f,a,b)
        calculate J2 and J4; J=J4 ; ERR=abs(J4-J2)

  2. Make a function DUMBINT that estimates the integral J=ab f(x)dx and the integration error by subdividing the interval (a,b) into N sub-intervals and using OPEN24 for each sub-interval:
    function [J,err] = DUMBINT (f,a,b,N)
       subdivide (a,b) into N sub-intervals (ai,bi)
       for each sub-interval calculate the sub-integral and sub-error [Si,ei] = OPEN24("f",ai,bi)
       calculate the global-integral J=∑iSi and global-error err=√(∑iei2)
    
  3. Using your program calculate the integral
    -100100  [1/1+x2] dx
    and estimate the number of points (integrand evaluations) necessary to achieve the accuracy of err = 0.0001.

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