The Simpson’s Method
1.Information About the Simpson’s
Method
The Newton-Cotes (given as follows)
formulas are generally unsuitable for use over large integration intervals
since high degree formulas would be required for the use over such intervals
and the values of the coefficients in these formulas are difficult to obtain.
Also the Newton Cotes formulas are based on interpolatory polynomials that use
equally spaced nodes, a procedure that is inaccurate over large intervals
because of the oscillatory nature of high degree polynomials.
n=1 Trapezoidal Rule
where 
n=2 Simpson’s Rule
where 
n=3 Simpson’s Three
Eighths Rule
where 
But here a piecewise approach in numerical analysis
that uses low order Newton-Cotes given above will be explained. This technique
is most often applied in practice.
To generalize this technique, choose an even
integer n. Subdivide the interval [a,b] into n subintervals and apply Simpson’s
Rule (provided above) on each consecutive pair of intervals.
With
and 
for each j = 0,1,……,n, we have
for some
with 
provided that 
. Using the fact that for each j=1,2,……,n/2-1, 
appears in the term
corresponding to the interval 
and also in the term
corresponding to the interval 
,this reduces to
where the last term is the error term. By the Intermediate Theorem this error term can be written as
where 
.
These derivations produce the following result
Let 
, n be even, 
and 
for each j = 0,1,……,n. There exists a 
for which the Composite
Simpson’s rule for n subintervals can be written with its error term as
The following algorithm uses the above expression to approximate the value of integral given on the left side of this expression.
INPUT
endpoints a,b; even positive integer n.
OUTPUT
approximation XI to the given integral.
Step 1: Set
Step 2: Set
XI0=f(a)+f(b)
XI1=0; (Summation of 
)
XI2=0. (Summation of 
)
Step 3:
For i=1,2,……,n-1 do Steps 4 and 5
Step 4:Set X=a+jh
Step 5: If i is even then set XI2=XI2+f(X)
else set XI1=XI1+f(X).
Step 6:
Set XI=h(XI0+2XI2+4XI1)/3.
Step 7:
OUTPUT(XI);
STOP.
This procedure can be applied to lower order formulas.
Composite Trapezoidal Rule is an example for this. You can see this if you
click on the Info Button when the Trapezoidal Rule is checked.
2.Adaptive Quadrature
The composite formulas require the use of equally spaced nodes.
For many problems this is not an important restriction, but it is an
inappropriate when integrating a function on an interval that contains both
regions with large functional variations and regions with small functional
variation. If the approximation error is to be evenly distributed, a smaller
step size is required for large variation regions than for those with less
variation. An efficient technique for this type problem can distinguish the
amount of functional variation and adapt the step size to the varying
requirements of the problem. These methods are known as Adaptive Quadratures
Methods. These methods can be applied to any composite rules with succes. In
our program we use these techniques for Trapezoidal and Simpson Rules and
Gaussian Quadratures.
The technique we use can be explained
as follows:
The sum of the approximation on the
sunintervals approximates the value of the integral on the whole integral. If the
approximation on one of the subintervals fails to be within the tolerance
epsilon/2, that subinterval is itself subdivided and its two subintervals analyzed
to determine if the approximation on each subinterval is accurate within
epsilon/4. This halving procedure is continued until each portion is within the
required tolerance. Although problems can be constructed for which this
tolerance will never be met but the this technique is successful for most of
the problems, because each subdivision typically increases the accuracy of the approximation
by a factor of 15 while requiring an increased accuracy factor of only 2.
3. Numerical Results and
Conclusions
f(x)=sin(x) on the interval [0,p]
The following figure is plot of this
function on the provided interval.
Figure.1: The plot of f(x)=sin(x)
on the provided interval.
The following table gives the results of the
integration of this function on the provided interval with different methods.
|
Trapezoidal |
Simpson’s |
Gaussian Quad. |
Romberg |
Epsilon |
|
1.9741510 |
2.0045765 |
1.9999906 |
1.9999999 |
0.01 |
|
1.9969667 |
2.0026127 |
1.9999999 |
1.9999999 |
0.001 |
|
1.9996646 |
2.0000007 |
1.9999999 |
1.9999999 |
0.0001 |
Table.1: The results for
f(x)=sin(x) on the interval [0, p]
f(x)=(100/x2)sin(10/x)
on the interval [1,3]
The following figure is the plot of this function on the provided
interval.
Figure.2: The plot of f(x)=(100/x2)sin(10/x)
on the interval [1,3]
The following table gives
the results of the integration of this function on the provided interval with
different methods.
|
Trapezoidal |
Simpson’s |
Gaussian Quad. |
Romberg |
Epsilon |
|
-1.3907865 |
-1.4294466 |
-1.4258159 |
-1.4260247 |
0.01 |
|
-1.4226971 |
-1.4258122 |
-1.4226024 |
-1.4260247 |
0.001 |
|
-1.4260349 |
-1.4226021 |
-1.4260241 |
-1.4260247 |
0.0001 |
Table.2: The results for
f(x)=(100/x2)sin(10/x) on the interval [1, 3]
f(x)=x2ex on the interval [0,1]
The following figure is plot of this
function on the provided interval.
Figure.3: The plot of f(x)=x2ex
on the interval [0,1]
The following table gives the results of the
integration of this function on the provided interval with different methods.
|
Trapezoidal |
Simpson’s |
Gaussian Quad. |
Romberg |
Epsilon |
|
0.1624884 |
0.1607224 |
0.1606027 |
0.1606280 |
0.01 |
|
0.1614894 |
0.1607224 |
0.1606027 |
0.1606280 |
0.001 |
|
0.1606369 |
0.1606103 |
0.1606027 |
0.1606279 |
0.0001 |
Table.3: The results for f(x)=x2ex
on the interval [0,1]
f(x)=cos(1/x) on the interval [0.1,2]
The following figure is plot of this
function on the provided interval.
Figure.4: The plot of
f(x)=cos(1/x) on the provided interval.
The following table gives the results of the
integration of this function on the provided interval with different methods.
|
Trapezoidal |
Simpson’s |
Gaussian Quad. |
Romberg |
Epsilon |
|
0.5814277 |
0.6739107 |
0.6528763 |
0.6738292 |
0.01 |
|
0.5814277 |
0.6737815 |
0.6738920 |
0.6738321 |
0.001 |
|
0.6736501 |
0.6738119 |
0.6738800 |
0.6738321 |
0.0001 |
Table.4: The results for
f(x)=cos(1/x) on the interval [0.1, 2]
f(x)=e2xsin(3x) on the interval [1,3]
The following figure is plot of this
function on the provided interval.
Figure.5: The plot of f(x)=e2xsin(3x)
on the interval [1,3]
The following table gives the results of the
integration of this function on the provided interval with different methods.
|
Trapezoidal |
Simpson’s |
Gaussian Quad. |
Romberg |
Epsilon |
|
108.5886632 |
108.5560839 |
108.5552994 |
108.5552994 |
0.01 |
|
108.5555266 |
108.5553434 |
108.5552994 |
108.5552994 |
0.001 |
|
108.5553235 |
108.5552971 |
108.5552812 |
108.5552812 |
0.0001 |
Table.5: The results for f(x)=e2xsin(3x)
on the interval [1,3]
By
looking at the above results, it can be said that all of the four methods are
good at calculating the integrals of oscillating functions. But the effect of
epsilon on different methods is different. The Trapezoidal and the Simpson’s
Rules are not giving sufficiently accurate results for epsilon bigger than
0.001. But the methods of Romberg Integration and Quassian Quadratures are
giving nearly same results and these results are better than the ones obtained
from the other two methods especially for bigger values of epsilon. The source
of this difference is the composite rules. The approximations done by
Trapezoidal and Simpson’s Methods for one subinterval are not as accurate as
the ones that are done by the other two methods. This also causes more levels
of halving for the Trapezoidal and Simpson’s Rules. Then although we didn’t
provide the exact numbers here, the time that is required for the calculation
of the integrals of the functions given above can be ordered as follows:
4. References
[1] Numerical Analysis, R.L. Burden, J.D.
Faires, PWS Publishing Company, Boston 1993.