Simpsons rule is a numerical integration method which makes use of parabolas. This provides a closer approximation to numerical integrals when compared to Left, Right Riemann sums and mid point rule which use rectangles and Trapezoidal rule which measures the area of trapezoids for approximation. This rule is given its name after the English mathematician Thomas Simpson who popularized this method of numerical integration.As it is done for other numerical integration methods, the interval [a,b] is divided into n equal partitions of width = Δx. But in order to apply this rule, n has to be even. The area in any two successive partition is approximated to the area under the parabola, which passes through the three end points which define the partition. In the diagram shown below the area in the partitions [xi-1 , xn+1] are approximated to the area of the region under the parabola shown which passes through the points P1, P2 and P3.

The Simpson's Rule for approximation of a definite integral is derived by using the general equation for one such parabola and evaluating the area under the parabola applying fundamental theorem of Calculus. Then the sum of all such areas is found using the pattern observed from the area of one sub interval and is known as Simpson's Rule for numerical integration.
Let us look how Simpson's rule is stated, its proof and also solve few example problems.

Simpson's rule used for approximating a definite integral is as follows:

$\int_{a}^{b}f(x)dx$ ≈ S = $\frac{\Delta x}{3}$$[f(x_{0})+4f(x_{1})+2f(x_{2})+4f(x_{3})+........+2f(x_{n-2})+4f(x_{n-1})+f(x_{n})]$

where 'n' is even and Δx = $\frac{b-a}{n}$.
The rule can be remembered observing pattern found in the sum as the function values are alternatively multiplied by 4 and 2 except the function values at the boundary points which are taken as such.

Simpson's rule approximation is more efficient when compared to mid point rule or Trapezoidal Rule. The Es error involved in using Simpson's rule is given by
E ≤ $\frac{K(b-a)^{5}}{180n^{4}}$ and K is so chosen that | f 4 (x) | ≤ K where f 4 (x) denotes the fourth derivative of f(x).

The interval [a, b] is divided into n sub intervals of equal width h = Δx = $\frac{b-a}{n}$ where n is an even integer.

On each consecutive pair of intervals the curve y = f(x) ≥ 0 is approximated to a parabola as shown in the diagram below. If we consider the general point on f(x) as Pi (xi, yi) corresponding to xi then yi = f(xi). We take a general parabola to pass through the points Pi, Pi+1 and Pi+2.



Even Partitions and Parabolas
Area under the first Parabola after the
horizontal shift.

Let us first find the equation of the parabola passing through the points P0, P1 and P2. We can assume x0 = -h, x1 = 0 and
x2 = h. The general form of the equation to a parabola is y = Ax2 + Bx + C. Thus applying the fundamental theorem of calculus, the area under the parabola from -h to h is,
$\int_{-h}^{h}(Ax^{2}+Bx+C)dx$ = 2$\int_{0}^{h}(Ax^{2}+C)dx$ Ax + C is an even function and Bx is an odd function
= 2$[A$$\frac{x^{3}}{3}$$+Cx]_{0}^{h}$

= 2$[A$$\frac{h^{3}}{3}$$+Ch]$

= $\frac{h}{3}$$[2Ah^{2}+6C]$

The parabola Ax2 + Bx + C = 0 contains the points (-h, y0), (0, y1) and (h, y2)
Thus we have y0 = Ah2 - Bh + C ...................(1)
y1 = C ....................(2)
y2 = Ah2 + Bh + C ...................(3)
Solving these three equations it can be shown
$2Ah^{2}+6C$ = y0 + 4y1 + y2.
Thus the area under the parabola passing through P0, P1 and P2 is

$\frac{h}{3}$ (y0 + 4y1 + y2)

As a horizontal shift of points P0, P1 and P2 resulting in x coordinates x0, x1 and x2 will not alter the area under the parabola passing through P0, P1 and P2 and it is still = $\frac{h}{3}$ (y0 + 4y1 + y2).
Similarly, the area under the parabola through P2, P3 and P4 from x = x2 to x = x4 is equal to
$\frac{h}{3}$ (y2 + 4y3 + y4)

When all the areas under all parabolas is computed using the above patter and summed up we have
$\int_{a}^{b}f(x)dx$ = $\frac{h}{3}$ (y0 + 4y1 + y2) + $\frac{h}{3}$ (y2 + 4y3 + y4) + $\frac{h}{3}$ (y4 + 4y5 + y6) + ...........

+ $\frac{h}{3}$ (yn-2 + 4yn-1 + yn)

= $\frac{h}{3}$ (y0 + 4y1 + 2y2 + 4y3 + 2y4 + ............ + 2yn-2 + 4yn-1 + yn)

= $\frac{b-a}{3n}$ (y0 + 4y1 + 2y2 + 4y3 + 2y4 + ............ + 2yn-2 + 4yn-1 + yn)

which is indeed known as Simpson's Rule.
Even though the approximation is derived for f(x) ≥ 0, it is a reasonable approximation whenever f(x) is continuous.

Solved Examples

Question 1: Use Simpson's rule to approximate the definite integral for the specified value of n.
    $\int_{0}^{4}\sqrt{x}sin(x)dx$      n =8.
Solution:
 
Substituting f(x) = $\sqrt{x}sin(x)$, n =8 and Δx = $\frac{4-0}{8}$ = 0.5 in Simpson's Rule,

    $\int_{0}^{4}\sqrt{x}sin(x)dx$ ≈ S8 = $\frac{\Delta x}{3}$ [f(0) + 4f(0.5) + 2f(1) + 4f(1.5) + 2f(2) + 4f(2.5) + 2f(3) + 4f(3.5)
                                                                                       + f(4)]
                                                         = $\frac{0.5}{3}$ [0 + 1.356 + 1.683 + 4.8868 + 2.5718 + 3.7852 + 0.4888   -2.6252 -1.5136]
                                                         = 0.1667 [ 10.6328] = 1.7721
 

Question 2: Approximate the definite integral using Simpson's rule for n = 4 and n = 10.

    $\int_{1}^{2}e^{\frac{1}{x}}dx$
Solution:
 
When n =4, Δx = $\frac{2-1}{4}$ = 0.25

    Simpson;s rule approximation is hence

    $\int_{1}^{2}e^{\frac{1}{x}}dx$  ≈ S4 = $\frac{\Delta x}{3}$ [f(1) + 4f(1.25) + 2f(1.5) + 4f(1.75) + f(2)]
                                                          = $\frac{0.25}{3}$ [2.7183 + 8.902 + 3.8954 + 7.0832 + 1.6487]
                                                          = 2.02063    (rounded to five decimal places)
   Simpson's approximation for n = 10
   Δx = 0.1
   $\int_{1}^{2}e^{\frac{1}{x}}dx$  ≈ S10 = $\frac{\Delta x}{3}$ [f(1) + 4f(1.1) +2f(1.2) +4f(1.3) + 2f(1.4) + 4f(1.5) + 2f(1.6)
                                                                                             + 4f(1.7) + 2f(1.8) + 4f(1.9) + f(2)]
                                                            = $\frac{0.1}{3}$ [2.7183 + 9.9284 + 4.602 + 8.6324 + 4.0854 + 7.7908 + 3.7364 +
                                                                            7.2032 + 3.4858 + 6.7708 + 1.6487]
                                                            = 2.02007
The approximation got for n= 10 is closer to the actual value of the integral = 2.0201