Simpson's Rule

If we are given odd number of tabular points,i.e. $n$ is even, then we can divide the given integral of integration in even number of sub-intervals $[x_{2k}, x_{2k+2}].$ Note that for each of these sub-intervals, we have the three tabular points $x_{2k},x_{2k+1}, x_{2k+2}$ and so the integrand is replaced with a quadratic interpolating polynomial. Thus using the formula (13.3.3), we get,

$$\int\limits^{x_{2k+2}}_{x_{2k}}f(x)dx = \frac{h}{3}\left[y_{2k} + 4 y_{2k+1} + y_{2k+2} \right].$$

In view of this, we have

$$\int\limits^{b}_{a}f(x)dx = \int\limits^{x_2}_{x_0}f(x)dx + \int\limits^{x_4}_{x_2}f(x)dx +\c... ...\limits^{x_{2k+2}}_{x_{2k}}f(x)dx +\cdots + \int\limits^{x_{n}}_{x_{n-2}}f(x)dx$$

$$= \frac{h}{3}\left[(y_{0} + 4 y_{1} + y_{2})+(y_{2} + 4 y_{3} + y_{4})+ \cdots+ (y_{n-2} + 4 y_{n-1} + y_{n}) \right]$$

$$= \frac{h}{3}\left[y_{0} + 4 y_{1} + 2 y_{2}+ 4 y_{3} + 2 y_{4}+ \cdots+ 2 y_{n-2} + 4 y_{n-1} + y_{n} \right],$$

which gives the second quadrature formula as follows:

$$\int\limits^{b}_{a}f(x)dx = \frac{h}{3}\left[(y_{0} +y_{n})+ 4\times(y_{1} + y_{3}+ \cdots + y_{2k+1}+ \cdots+ y_{n-1})\right.$$

$$+ \left. 2\times (y_{2}+ y_{4}+ \cdots+ y_{2k}+\cdots+ y_{n-2})\right]$$

$$= \frac{h}{3}\left[(y_{0} +y_{n})+ 4\times\left(\sum_{i=1,\; i - od... ...1}y_{i}\right) + 2\times\left(\sum_{i=2,\; i - even }^{n-2}y_{i}\right)\right].$$ (13.3.5)

This is known as SIMPSON'S RULE.

Remark 13.3.3   An estimate for the error $E_2$ in numerical integration using the Simpson's rule is given by

$$E_2= - \frac{b-a}{180}\overline{\Delta^4 y},$$ (13.3.6)

where $\overline{\Delta^4 y}$ is the average value of the forth forward differences.

EXAMPLE 13.3.4   Using the table for the values of $y = e^{x^2}$ as is given in Example 13.3.2, compute the integral $\int\limits^{1}_{0}e^{x^2}dx,$ by Simpson's rule. Also estimate the error in its calculation and compare it with the error using Trapezoidal rule.
Solution: Here, $h=0.1, \; n=10,$ thus we have odd number of nodal points. Further,

$$y_0+y_{10}= 1.0+2.71828= 3.71828, \;\;\; \sum_{i=1,\; i - odd }^{9}y_{i}= y_1 +y_3+y_5+y_7+y_9=7.26845,$$

and

$$\sum_{i=2,\; i - even }^{8}y_{i}=y_2+y_4+y_6+y_8 =5.54412.$$

Thus,

$$\int\limits^{1}_{0}e^{x^2}dx= \frac{0.1}{3}\times[3.71828 + 4\times7.268361+ 2\times5.54412 ]=1.46267733$$

To find the error estimates, we consider the forward difference table, which is given below:

$x_i$	$y_i$	$\Delta y_i$	$\Delta^2 y_i$	$\Delta ^3 y_i$	$\Delta ^4 y_i$
0.0	1.00000	0.01005	0.02071	0.00189	0.00149
0.1	1.01005	0.03076	0.02260	0.00338	0.00171
0.2	1.04081	0.05336	0.02598	0.00519	0.00243
0.3	1.09417	0.07934	0.03117	0.00762	0.00320
0.4	1.17351	0.11051	0.3879	0.01090	0.00459
0.5	1.28402	0.14930	0.04969	0.01549	0.00658
0.6	1.43332	0.19899	0.06518	0.02207	0.00964
0.7	1.63231	0.26417	0.08725	0.03171
0.8	1.89648	0.35142	0.11896
0.9	2.24790	0.47038
1.0	2.71828

Thus, error due to Trapezoidal rule is,

$$E_1 = - \frac{1-0}{12}\overline{\Delta^2 y}$$

$$= - \frac{1}{12}\times\frac{0.02071+0.02260+0.02598+0.03117+0.03879+0.04969+0.06518+0.08725+0.11896}{9}$$

$$= -0.004260463.$$

Similarly, error due to Simpson's rule is,

$$E_2 = - \frac{1-0}{180}\overline{\Delta^4 y}$$

$$= - \frac{1}{180}\times\frac{0.00149+0.00171+0.00243+0.00328+0.00459+0.00658+0.00964}{7}$$

$$= -2.35873\times 10^{-5}.$$

It shows that the error in numerical integration is much less by using Simpson's rule.

EXAMPLE 13.3.5   Compute the integral $\int\limits^{1}_{0.05}f(x)dx$ , where the table for the values of $y=f(x)$ is given below:

$x$	0.05	0.1	0.15	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1.0
$y$	0.0785	0.1564	0.2334	0.3090	0.4540	0.5878	0.7071	0.8090	0.8910	0.9511	0.9877	1.0000

Solution: Note that here the points are not given to be equidistant, so as such we can not use any of the above two formulae. However, we notice that the tabular points $0.05, 0.10, 0,15$ and $0.20$ are equidistant and so are the tabular points $0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9$ and $1.0$ . Now we can divide the interval in two subinterval: $[0.05, 0.2]$ and $[0.2, 1.0]$ ; thus,

$$\int\limits^{1}_{0.05}f(x)dx=\int\limits^{0.2}_{0.05}f(x)dx+\int\limits^{1}_{0.2}f(x)dx$$

. The integrals then can be evaluated in each interval. We observe that the second set has odd number of points. Thus, the first integral is evaluated by using Trapezoidal rule and the second one by Simpson's rule (of course, one could have used Trapezoidal rule in both the subintervals).
For the first integral $h=0.05$ and for the second one $h=0.1$ . Thus,

$$\int\limits^{0.2}_{0.05}f(x)dx = 0.05\times\left[\frac{0.0785+0.3090}{2}+0.1564+0.2334\right]= 0.0291775,$$

$${\mbox{ and }} \;\; \int\limits^{1.0}_{0.2}f(x)dx = \frac{0.1}{3}\times \biggl[(0.3090+1.0000)+4\times(0.4540+0.7071+0.8910+ 0.9877) \biggr.$$

$$\hspace{1in} \biggl. +2\times( 0.5878+0.8090+ 0.9511)\biggr]$$

$$= 0.6054667,$$

which gives,

$$\int\limits^{1}_{0.05}f(x)dx= 0.0291775+0.6054667 = 0.6346442$$

It may be mentioned here that in the above integral, $f(x)=sin(\pi x/2)$ and that the value of the integral is $0.6346526$ . It will be interesting for the reader to compute the two integrals using Trapezoidal rule and compare the values.

EXERCISE 13.3.6  

1. Using Trapezoidal rule, compute the integral $\int\limits^{b}_{a}f(x)dx,$ where the table for the values of $y=f(x)$ is given below. Also find an error estimate for the computed value.

   1. $x$	a=1	2	3	4	5	6	7	8	9	b=10
      $y$	0.09531	0.18232	0.26236	0.33647	0.40546	0.47000	0.53063	0.58779	0.64185	0.69314

   2. $x$	a=1.50	1.55	1.60	1.65	1.70	1.75	b=1.80
      $y$	0.40546	0.43825	0.47000	0.5077	0.53063	0.55962	0.58779

   3. $x$	a = 1.0	1.5	2.0	2.5	3.0	b = 3.5
      $y$	1.1752	2.1293	3.6269	6.0502	10.0179	16.5426

2. Using Simpson's rule, compute the integral $\int\limits_a^b f(x) dx.$ Also get an error estimate of the computed integral.
   1. Use the table given in Exercise 13.3.6.1b.

   2. $x$	a = 0.5	1.0	1.5	2.0	2.5	3.0	b = 3.5
      $y$	0.493	0.946	1.325	1.605	1.778	1.849	1.833

3. Compute the integral $\int\limits^{1.5}_{0}f(x)dx$ , where the table for the values of $y=f(x)$ is given below:
   

   $x$	0.0	0.5	0.7	0.9	1.1	1.2	1.3	1.4	1.5
   $y$	0.00	0.39	0.77	1.27	1.90	2.26	2.65	3.07	3.53