Subject : The following demonstrates the convergence of the Trapezoidal method of estimating integrals of continuous functions.
Authors : Nathan Collier, Autar Kaw, Loubna Guennoun
Date : 27 October 2005

Introduction

Simpson's rule is based on the Newton-Cotes formula that if one approximates the integrand of the integral by an nth order polynomial, then the integral of the function is approximated by the integral of that nth order polynomial. Integration of polynomials is simple and is based on the calculus.  Simpson's 1/3rd rule is the area under the curve where the function is approximated by a second order polynomial.   [click here for textbook notes] [click here for power point presentation].

The following simulation illustrates the convergence of Simpson's 1/3rd rule of integration. This section is the only section where the user interacts with the program.  The user enters any function f(x), the lower and upper limit of the integration, and the maximum number of segments, n. The program will display the true error, the absolute relative percentage true error, the approximate error, the absolute relative percentage approximate error, and the least number of significant digits correct in the approximation.

Inputs

Integrand f[x] = 0

f[x_] := 300.0*x/(1.0 + Exp[x]);

Lower and upper limit of the integral, a and b respectively.

a = 0.0;
b = 10.0;

curve = 
        Plot[f[x], {x, a, b}, 
          PlotLabel -> "Entered function on given interval", 
          TextStyle -> {FontSize -> 
                11}];
General::spelll  Off;

Maximum number of segments. Note that this number must be even.

(Subscript[n, maximum]) = 40;

True Solution

This is the solution found by Mathematica

Actual = Integrate[f[x], {x, a, b}]

Value of integral as a function of segments

The following procedure determines the value of the integral from a to b using Simpson's 1/3rd rule with n segments.

Array[AV, (Subscript[n, maximum])];

For[i = 2, i ≤ (Subscript[n, maximum]), 
      i = i + 2, 
h = (b - a)/
          i; 
sum1 = 
        Sum[f[a + jj*h], {jj, 1, i - 1, 2}]; 
sum2 = 
        Sum[f[a + jj*h], {jj, 2, i - 2, 2}]; 
AV[i] = 
        h/3*(f[a] + 4*sum1 + 2*sum2 + f[b])]

Absolute true error

Array[Et, (Subscript[n, maximum])];

For[i = 2, i <= (Subscript[n, maximum]), i = i + 2, 
      Et[i] = Abs[Actual - AV[i]]]

Absolute relative true error

Array[ϵt, (Subscript[n, maximum])];

For[i = 2, i <= (Subscript[n, maximum]), 
      i = i + 2, ϵt[i] = Abs[Et[i]/Actual*100]]

Absolute approximate error

Array[Ea, (Subscript[n, maximum])];

For[i = 2, i <= (Subscript[n, maximum]), i = i + 2, 
      If[i ≤ 2,  Ea[i] =  0, 
        Ea[i] = Abs[AV[i] - AV[i - 2]]]]

Absolute relative approximate error

Array[ϵa, (Subscript[n, maximum])];

For[i = 2, i <= (Subscript[n, maximum]), i = i + 2, 
      If[i ≤ 2, ϵa[i] = 0, ϵa[i] = 
          Abs[Ea[i]/AV[i]*100]]]

Significant digits at least correct

Array[sigdig, (Subscript[n, maximum])];

For[i = 2, i <= (Subscript[n, maximum]), i = i + 2, 
      If[(ϵa[i] ≥ 5)  ||  (i ≤ 
              2), sigdig[i] = 0, 
        sigdig[i] = 
          Floor[(2 - Log[10, Abs[ϵa[i]/0.5]])]]]

Graphs

xrplot = Table[AV[i], {i, 2, (Subscript[n, maximum]), 2}];

ListPlot[xrplot,PlotJoined -> True, PlotRange -> All, AxesOrigin -> {1, Min[xrplot]}, PlotLabel -> "Estimated root as a function of number \
of iterations"];

Etplot = Table[Et[i], {i, 2, (Subscript[n, maximum]), 2}];

ListPlot[Etplot,PlotJoined -> True, PlotRange -> All, AxesOrigin -> {1, Min[Etplot]}, PlotLabel -> "Absolute true error as a function of \
number of iterations"];

ϵtplot = 
        Table[ϵt[i], {i, 2, (Subscript[n, maximum]), 2}];

ListPlot[ϵtplot,PlotJoined -> True, PlotRange -> All, AxesOrigin -> {1, Min[ϵtplot]}, PlotLabel -> "Absolute relative true error as a \
function of number of iterations"];

Eaplot = Table[Ea[i], {i, 2, (Subscript[n, maximum]), 2}];

ListPlot[Eaplot,PlotJoined -> True, PlotRange -> All,AxesOrigin -> {1, Min[Eaplot]}, PlotLabel -> "Absolute approximate error as a function \
of number of iterations"];

ϵaplot = 
        Table[ϵa[i], {i, 2, (Subscript[n, maximum]), 2}];

ListPlot[ϵaplot,PlotJoined -> True, PlotRange -> All, AxesOrigin -> {1, Min[ϵaplot]}, PlotLabel -> "Absolute relative approximate error as a \
function of number of iterations"];

sigdigplot = Table[sigdig[i], {i, 2, (Subscript[n, maximum]), 2}];

<< Graphics`Graphics`

BarChart[sigdigplot];