Approximate Integration

Sometimes it is necessary to obtain a numerical estimate of the definite integral of a function f(x) between x=a and x=b (where b>a). This integral is represented by the area bounded by the curve y=f(x), the x-axis and the ordinates at x=a and x=b. One reasonably accurate method of estimating this area, and hence this integral, is to use Simpson's Rule. Here the area is approximated by

(h/3)(y₀ + 4y₁ + y₂) + (h/3)(y₂ + 4y₃ + y₄) + ... + (h/3)(y_2n-2 + 4y_2n-1 + y_2n)

The area is divided into 2n strips of width h=(b-a)/2n and the ordinates are given by y₀=f(a), y₁=f(a+h), y₂=f(a+2h) etc. The method uses a quadratic of the form AX²+BX+C to approximate the function between three adjacent points on the curve (i.e. between two adjacent strips).