Formulae With Equally Spaced Points

The formula in the above equation can be written in a compact form by defining $h = \frac{ b - a}{N}$. What this means is the interval $[a,b]$ is divided into N equal parts of length $h$ and the integral in each part is approximated by $h$ times the value of the function at the midpoint ( the first equation ) or the average of the function at two ends of the part ( the second equation ). These formulae are often called as extended midpoint rule and extended trapizoidal rule respectively. The meaning of these terms is obvious. The formula is exact in the limit $N \rightarrow \infty$. On a computer, however, we cannot take the limit $N \rightarrow \infty$ and we have to estimate the integral for some reasonably large value of N. This is depicted in Figure below. In this figure the trapizoidal and midpoint rules are shown by dashed and dotted lines respectively.

Figure 2: Geometric definition of Integral

In fact, there is an exact expression for the integral in terms of the sum given above plus correction terms which depend on the higher derivatives of the integrand. :

$$I = \frac{b-a}{N} \Big [ \frac{1}{2} f(a) + \sum_{i = 1, N-1} f \Big ... ...2} f(b) \Big ] - \sum_i \frac{B_{2i} h^{2i}}{ (2i)! } ( f^{2i}(b) - f^{2i}(a) )$$

$$= I^t_N - \sum_i \frac{B_{2i} h^{2i}}{ (2i)! } ( f^{2i}(b) - f^{2i}(a) )$$ (14)

where $B_i$ are the Bernaulli numbers. This is the Euler-Maclaurin summation formula for the integration. So, in principle, if one knows all the deriatives of the function, one has an analytic expression for the integral. Usually, this is not the case so one has to limit the summation to finite number of derivatives. In fact, most of the time we may not know any derivatives of the integrand at the end points.

Clearly, the correction terms go to zero ( provided all the higher derivatives of the function are finite ) in the limit $h = \frac{(b-a)}{N}$ going to zero. We can take advantage of the formula by computing the first term for two or more values of N and use the result to eliminate some of the leading terms of the correction series. In particular, if we compute the integral for N and 2N intervals then it is straight forward to show that

$$I = \frac{4}{3} I^t_{2N} - \frac{1}{3} I^t_{N} + {\cal O} (\frac{1}{N^4} )$$ (15)

So, with the two computations we have got a new formula which is expected to be much more accurate. This is the extended Simpson's rule. We can write it in a more familiar form by defining $h = \frac{b-a}{2N}$,

$$I^S_{2N} = h \Big [ \frac{1}{3} f(a) + \frac{4}{3} f(a+h) + \frac{2}{3} f(a+2h) \cdots \frac{1}{3} f(b)$$ (16)

The trapizoidal and Simpson rules are called closed formulae because the end-points are included in the formula. The midpoint rule is an open formula. Sometimes open formula is useful if the integrand has an integrable singularity at the end-point(s). The closed formula does not work in such cases because the function is singular at the end point(s).

We can obtain a formula similar ( in spirit ) to the Simpson formula from the midpoint rule. The trick is to triple the number of points here since by doubling we cannot get advantage of previous computations. We then have

$$I = \frac{9}{8} I^m_{3N} - \frac{1}{9} I^m_N + {\cal O} (\frac{1}{N^4} )$$ (17)

where $I^m_N$ is the estimate of the integral using extended midpoint rule;

$$I = \frac{b-a}{N} \sum_{i = 1,N} f\Big ( a + \frac{(i - {1 \over 2} )(b-a)}{N} \Big )$$ (18)

• Program 1: Write down a subroutine which will compute an integral $\int_a^b dx f(x)$ using trapizoidal rule. The input should be the lower and upper limits ( a and b ), the function f which is defined as an external function ( to be defined separately ) and output should be the estimate of the integral and an error estimate. The procedure adopted should be that begining with number of intervals equal to 1, the program should double the number of intervals every time and compute the integral. You can choose maximum number of intervals to be $2^N$ with $N$ some integer ( say between 5 and 10 ). Error extimate will then be the difference between the integral for $2^{N-1}$ and $2^N$ intervals. Other possibility is to keep on increasing the number of intervals till some desired accuracy is reached ( i.e. $\vert ( I_{2^{N-1}} - I_{2^N} ) / I_{2^N} \vert < \epsilon$ ).
• Program 2: Write the subroutine for computing the integral using Simpson's rule. Note that you can modify the program developed in the previous problem to do Simpson's rule. But then N should be greater than 2 ( at least ). You will then have to use the result for $2^{N-1}$ and $2^N$ intervals to compute the integral using Simpson's rule. You can also have error extimate as in the previous problem.
• Problem 3: One can use extrapolation technique to get a better estimate of the integral using trapizoidal rule for $2^N$ intervals with different values of N. From Euler-Maclaurin series we know that the correction terms form a series in powers of $1/N^2$. Considering the integral estimated from trapizoidal rule as a function of the variable $1/N^2$, we can estimate the integral by extrapolation to $1/N^2 = 0$. One can use polynomeal interpolation formula.

There are a number of other integration formulae ( both open and closed type ) using intervals of equal spacing. We shall not be considering them as they are more complicated to program and they do not seem to have much advantage over trapizoidal, midpoint or Simpson rules described above.

Given an integral to be done on computer, one is not sure about the number of points to be used in the calculation. In that case, one can use certain number of integration points and check the accuracy by doubling the number of points. Note that, if one uses trapizoidal rule, half of the points at which the function is needed for second calculation are already computed in the first step. Using the two results, we have a refined estimate of the answer ( which is actually the Simpson rule ). One should have the difference between the second computation and the refined estimate to be tolerably small. If it is not, one needs to double the number of points again and repeat the procedure. One can have an integration subroutine using this algorithm.