## Simpson's Rule

Simpson's rule approximates the integral of a function*f(x)*on the closed and bounded interval

*[a, a+h]*of length

*h > 0*by the integral on

*[a, a+h]*of the quadratic passing through the points

*(a, f(a))*,

*(a+h/2, f(a+h/2))*and

*(a+h, f(a+h))*. The composite Simpson's rule is used to approximate the integral of a function

*f(x)*over a closed and bounded interval

*[a, b]*where

*a < b*, by decomposing the interval

*[a, b]*into

*n > 1*subintervals of equal length

*h = (b - a) / n*, then adding the results of applying the Simpson's rule to each subinterval. By abuse of language both the composite Simpson's rule and Simpson's rule sometimes are referred to simply as Simpson's rule.

Let

*∫*be the integral of

_{a}^{b}f( x ) dx*f(x)*over the closed and bounded interval

*[a,b]*, and let

*S*be the result of applying the Simpson's rule with

_{h}(f)*n*subintervals of length

*h*, i.e.

*S*=

_{h}(f)*(h/6) [ f(a) + 4f(a+h/2) + 2f(a+h) + ··· + 2f(b-h) + 4f(b-h/2) + f(b) ]*.

An immediate consequence of the Euler-Maclaurin summation formula relates

*∫*and

_{a}^{b}f( x ) dx*S*

_{h}(f)*S*=

_{h}(f)*∫*+

_{a}^{b}f( x ) dx*(h*

+ ··· + K h,

^{4}/2880) [ f'''(b) - f'''(a) ] - (h^{6}/96768) [ f^{(5)}(b) - f^{(5)}(a) ]+ ··· + K h

^{ 2p-2}[f^{(2p-3)}(b) - f^{(2p-3)}(a) ] + O(h^{ 2p})where

*f'''*,

*f*, and

^{(5)}*f*are the third, fifth and (p-3)rd derivatives of

^{(2p-3)}*f*and

*K*is a constant.

The last term,

*O(h*is important. Given an infinitely differentiable function in which the first

^{ 2p})*2p-3*derivatives vanish at both endpoints of the interval of integration, it is not true that

*S*=

_{h}(f)*∫*but rather what the theorem says is that

_{a}^{b}f( x ) dx*lim*-

_{h→0}| ( S_{h}(f)*∫*) /

_{a}^{b}f( x ) dx*h*<

^{2p}|*M*,

*M > 0*.

If

*f*is at least four times differentiable on the interval

*[a,b]*, then applying the mean-value theorem to

*S*-

_{h}(f)*∫*=

_{a}^{b}f( x ) dx*(h*

+ ··· + K h

^{4}/2880) [ f'''(b) - f'''(a) ] - (h^{6}/96768) [ f^{(5)}(b) - f^{(5)}(a) ]+ ··· + K h

^{ 2p-2}[f^{(2p-3)}(b) - f^{(2p-3)}(a) ] + O(h^{ 2p})yields the standard truncation error expression

*S*-

_{h}(f)*∫*=

_{a}^{b}f( x ) dx*(h*, for some point

^{4}/2880) (b-a) f^{(4)}(c)*c*where

*a ≤ c ≤ b*.

A corollary of which is that if

*f*for all

^{(4)}(x) = 0*x*in

*[a,b]*, i.e. if

*f(x)*is a cubic, then Simpson's rule is exact.

The Euler-Maclaurin summation formula also shows that usually

*n*should be chosen large enough so that

*h = (b - a) / n < 1*. For example, if

*h = 0.1*then

*S*=

_{0.1}(f)*∫*+

_{a}^{b}f( x ) dx*3.5·10*

^{-8}[ f'''(b) - f'''(a) ] - 1.033·10^{-11}[ f^{(5)}(b) - f^{(5)}(a) ] + ···*h = 0.01*then

*S*=

_{0.01}(f)*∫*+

_{a}^{b}f( x ) dx*3.5·10*

^{-12}[ f'''(b) - f'''(a) ] - 1.033·10^{-17}[ f^{(5)}(b) - f^{(5)}(a) ] + ···*h = 10*then

*S*=

_{10}(f)*∫*+

_{a}^{b}f( x ) dx*3.47 [ f'''(b) - f'''(a) ] - 10.33 [ f*

^{(5)}(b) - f^{(5)}(a) ] +···*f(x)*is a cubic, then

*n*may be chosen to be 1.

Besides the truncation error described above, the algorithm is also subject to the usual round-off errors and errors due to number of significant bits in the machine representation of floating point numbers.

The source code below is programmed in double precision, but the number of significant bits can easily be extended to extended precision by changing

*double*to

*long double*and by affixing an

*L*to any floating point number e.g.

*0.5*to

*0.5L*.

Mathematically addition is associative,

*(a+b)+c = a + (b+c)*, but if intermediate results are rounded to machine precision, machine addition is not necessarily associative. E.g. consider a decimal machine with 3 significant digits, then

( ( 0.400 + 0.400 ) + 0.400 ) + 122 = 123 > 122 = 0.400 + ( 0.400 + ( 0.400 + 122 ) ).

In order to reduce the effect of intermediate result round-off errors when adding a series of floating point numbers, there are two versions of the rectangular rule, one which adds left to right
*S*=

_{h}(f)*h/6 [ (((···( ( f(a) + 4f(a+h/2) ) + 2f(a+h) ) + ··· + 4f(b-h/2) ) + f(b) ]*

*S*=

_{h}(f)*h/6 [ f(a) + ( 4f(a+h/2) + ( 2f(a+h) + ··· + (4f(b-h/2) + f(b) )···))) ]*.

### Function List

- double Simpsons_Rule_Sum_LR( double a, double h, int n, double (*f)(double) )

Integrate the user supplied function (*f)(x) from a to a + nh where a is the lower limit of integration, h > 0 is the length of each subinterval, and n > 0 is the number of subintervals. The sum is performed from left to right.

- double Simpsons_Rule_Sum_RL( double a, double h, int n, double (*f)(double) )

Integrate the user supplied function (*f)(x) from a to a + nh where a is the lower limit of integration, h > 0 is the length of each subinterval, and n > 0 is the number of subintervals. The sum is performed from right to left.

- double Simpsons_Rule_Tab_Sum_LR( double h, int n, double f[ ] )

Integrate the function f[ ] given as an array of dimension n whose i^{th}element is the function evaluated at the midpoint of the i^{th}subinterval, where h > 0 is the length of each subinterval, and n > 0 is the number of subintervals. The sum is performed from left to right.

- double Simpsons_Rule_Tab_Sum_RL( double h, int n, double f[ ] )

Integrate the function f[ ] given as an array of dimension n whose i^{th}element is the function evaluated at the midpoint of the i^{th}subinterval, where h > 0 is the length of each subinterval h, and n > 0 is the number of subintervals. The sum is performed from right to left.

*C* Source Code

- The file, simpsons_rule.c, contains the versions of Simpsons_Rule_Sum_LR( ) and Simpsons_Rule_Sum_RL( ) written in
*C*.

- The file, simpsons_rule_tab.c, contains the versions of Simpsons_Rule_Tab_Sum_LR( ) and Simpsons_Rule_Tab_Sum_RL( ) written in
*C*.