Numerical integration with Simpson's rule

Initialization cell, containing programming

<< Graphics`FilledPlot`
<< Graphics`Colors`

Left endpoint approximation

\(\(\(
\)\(\(leftendpoint[n_]\ \  := \ 
      Module[\ {X, \ dX}, 
X\  = \ 
          Table[a\  + \ k*\((b - a)\)/n, \ {k, \ 0, \ 
              n}]; 
Y\  = \ \ \ Map[f, \ 
            X]; 
verticals\  = \ 
          Table[\ Line[\ {
{X[\([k]\)], \ 
                  0}, \ {X[\([k]\)], \ 
                  Y[\([k]\)]}
}]\ \ , \ \ {k, \ 1, \ 
              n + 1}]; 
plot1\  = \ 
          Plot[f[x], \ {x, \ a, \ b}, \ 
            PlotStyle\  -> {AbsoluteThickness[2], \ Blue}, 
            PlotRange -> \ All, \ \ DisplayFunction\  -> \ 
              Identity]; 
functiontable\  = \ \ Table[\ {\ 
              Y[\([k]\)], \ \ \ \ X[\([k]\)]\  <= \ x\  && \ \ x\  <= \ 
                  X[\([k + 1]\)]}, \ \ {k, \ 1, \ n}\ ]; 
g[
            x_]\  := \ \ Piecewise[
            functiontable]; 
plot2\  = \ 
          FilledPlot[
            g[x], \ \ {x, \ a, \ 
              b}, 
Fills -> \ {GrayLevel[ .8]}, 
            PlotRange -> \ All, \ 
\ 
            DisplayFunction\  -> \ 
              Identity]; 
plot3\  = \ 
          Graphics[verticals]; 
leftsum\  = \ 
          N[Sum[Y[\([k]\)], \ \ {k, \ 1, \ n}]\ *\((b - a)\)/
                n]; 
label\  = 
          leftsum\ \ \ \ "  =   left endpoint sum  "\ \ ; \ 
Show[\ plot2, \ \ plot3, 
          plot1, 
PlotLabel\  -> \ 
            label, \ 
\ 
          DisplayFunction\  -> \ $DisplayFunction]
]\)\(
\)\(
\)\(
\)\(
\)
    \)\)\)

Midpoint approximation

\(\(\(
\)\(midpoint[n_]\ \  := \ 
      Module[\ {X, \ Y, \ Z, \ midpoints, verticals, \ plot1, plot2, \ 
          plot3, \ functiontable, \ g, \ }, 
X\  = \ 
          Table[a\  + \ k*\((b - a)\)/n, \ {k, \ 0, \ 
              n}]; 
midpoints\  = \ 
          Table[\((X[\([k]\)]\  + \ X[\([k + 1]\)])\)/2, \ \ {k, \ 1, \ 
              n}]; 
Y\  = \ \ \ Map[f, \ 
            midpoints]; 
Z\  = \ 
          Map[f, \ X]; 
verticals\  = \ 
          Table[\ Line[\ {
{X[\([k]\)], \ 
                  0}, \ {X[\([k]\)], \ 
                  Z[\([k]\)]}
}]\ \ , \ \ {k, \ 1, \ 
              n + 1}]; 
plot1\  = \ 
          Plot[f[x], \ {x, \ a, \ b}, \ 
            PlotStyle\  -> {AbsoluteThickness[2], \ 
                Blue}, 
PlotRange -> \ 
              All, \ \ DisplayFunction\  -> \ 
              Identity]; 
functiontable\  = \ \ Table[\ {\ 
              Y[\([k]\)], \ \ \ \ X[\([k]\)]\  <= \ x\  && \ \ x\  <= \ 
                  X[\([k + 1]\)]}, \ \ {k, \ 1, \ n}\ ]; 
g[
            x_]\  := \ \ Piecewise[
            functiontable]; 
plot2\  = \ 
          FilledPlot[g[x], \ \ {x, \ a, \ b}, 
            Fills -> \ {GrayLevel[ .8]}, \ \ PlotRange -> \ All, 
            DisplayFunction\  -> \ 
              Identity]; 
plot3\  = \ 
          Graphics[verticals]; 
midsum\  = \ 
          N[Sum[Y[\([k]\)], \ \ {k, \ 1, \ n}]\ *\((b - a)\)/
                n]; 
label\  = 
          midsum\ \ \ \ "  =   midpoint sum  "\ \ ; \ 
Show[\ plot2, \ \ plot3, \ 
          plot1, 
PlotLabel\  -> \ 
            label, \ 
\ 
          DisplayFunction\  -> \ $DisplayFunction]
]\)\)\)

Right endpoint approximation

\(\(\(
\)\(rightendpoint[n_]\ \  := \ 
      Module[\ {X, \ Y, \ verticals, \ functiontable, \ g, \ plot1, \ 
          plot2, \ plot3, \ leftsum, \ label}, 
X\  = \ 
          Table[a\  + \ k*\((b - a)\)/n, \ {k, \ 0, \ 
              n}]; 
Y\  = \ \ \ Map[f, \ 
            X]; 
verticals\  = \ 
          Table[\ Line[\ {
{X[\([k]\)], \ 
                  0}, \ {X[\([k]\)], \ 
                  Y[\([k]\)]}
}]\ \ , \ \ {k, \ 1, \ 
              n + 1}]; 
plot1\  = \ 
          Plot[f[x], \ {x, \ a, \ b}, \ 
            PlotStyle\  -> {AbsoluteThickness[2], \ Blue}, 
            PlotRange -> \ All, \ \ \ DisplayFunction\  -> \ 
              Identity]; 
functiontable\  = \ \ Table[\ {\ 
              Y[\([k + 1]\)], \ \ \ \ X[\([k]\)]\  <= \ x\  && \ \ x\  <= \ 
                  X[\([k + 1]\)]}, \ \ {k, \ 1, \ n}\ ]; 
g[
            x_]\  := \ \ Piecewise[
            functiontable]; 
plot2\  = \ 
          FilledPlot[g[x], \ \ {x, \ a, \ b}, 
            Fills -> \ {GrayLevel[ .8]}, \ \ PlotRange -> \ All, 
            DisplayFunction\  -> \ 
              Identity]; 
plot3\  = \ 
          Graphics[verticals]; 
rightsum\  = \ 
          N[Sum[Y[\([k + 1]\)], \ \ {k, \ 1, \ n}]\ *\((b - a)\)/
                n]; 
label\  = 
          rightsum\ \ \ \ "  =   right endpoint sum  "\ \ ; \ 
Show[\ plot2, \ \ plot3, \ 
          plot1, 
PlotLabel\  -> \ 
            label, \ 
\ 
          DisplayFunction\  -> \ $DisplayFunction]
]\)\)\)

Trapezoid approximation

\(\(\(
\)\(trapezoid[n_]\ \  := \ 
      Module[\ {X, \ Y, \ \ midpoints, verticals, \ plot1, plot2, \ plot3, \ 
          functiontable, \ g, \ \ dX\ }, 
X\  = \ 
          Table[a\  + \ k*\((b - a)\)/n, \ {k, \ 0, \ 
              n}]; 
Y\  = \ \ \ Map[f, 
            X]; 
dX\  = \ \((b - a)\)/
            n; 
verticals\  = \ 
          Table[\ Line[\ {
{X[\([k]\)], \ 
                  0}, \ {X[\([k]\)], \ 
                  Y[\([k]\)]}
}]\ \ , \ \ {k, \ 1, \ 
              n + 1}]; 
plot1\  = \ 
          Plot[f[x], \ {x, \ a, \ b}, \ 
            PlotStyle\  -> {AbsoluteThickness[2], \ 
                Blue}, 
PlotRange -> \ 
              All, \ \ DisplayFunction\  -> \ 
              Identity]; 
functiontable\  = \ \ Table[\ {\ 
              Y[\([k]\)]*\ \((X[\([k + 1]\)] - x)\)\ /
                    dX + \ \ Y[\([k + 1]\)]*\((x\  - \ X[\([k]\)])\)/
                    dX\ , \ \ \ \ X[\([k]\)]\  <= \ x\  && \ \ x\  <= \ 
                  X[\([k + 1]\)]}, \ \ {k, \ 1, \ n}\ ]; 
g[
            x_]\  := \ \ Piecewise[
            functiontable]; 
plot2\  = \ 
          FilledPlot[g[x], \ \ {x, \ a, \ b}, 
            Fills -> \ {GrayLevel[ .8]}, \ 
PlotRange -> \ All, \ \ DisplayFunction\  -> \ 
              Identity]; 
plot3\  = \ 
          Graphics[
            verticals]; 

trapezoidsum\  = \ \(( .5*Y[\([1]\)]\  + \  .5*Y[\([n + 1]\)]\  + \ 
                Sum[Y[\([i]\)], \ {i, \ 2, \ n}])\)\ *\ 
            dX; 
label\  = 
          trapezoidsum\ \ \ \ "  =  trapezoid sum  "\ \ ; \ 
Show[\ plot2, \ \ plot3, \ 
          plot1, 
PlotLabel\  -> \ 
            label, \ 
\ 
          DisplayFunction\  -> \ $DisplayFunction]
]\)\)\)

Simpson approximation

\(\(\(
\)\(simpson[n_]\ \  := \ 
      Module[\ {X, \ Y, \ verticals, \ plot1, plot2, \ plot3, \ 
            functiontable, \ g, \ \ dX\ , \ 
            simpsonsum}, 
X\  = \ 
            Table[a\  + \ k*\((b - a)\)/n, \ {k, \ 0, \ 
                n}]; 
Y\  = \ \ \ Map[f, 
              X]; 
dX\  = \ \((b - a)\)/
              n; 
verticals\  = \ 
            Table[\ Line[\ {
{X[\([k]\)], \ 
                    0}, \ {X[\([k]\)], \ 
                    Y[\([k]\)]}
}]\ \ , \ \ {k, \ 1, \ 
                n + 1}]; 
plot1\  = \ 
            Plot[f[x], \ {x, \ a, \ b}, \ 
              PlotStyle\  -> {AbsoluteThickness[2], \ 
                  Blue}, 
PlotRange -> \ 
                All, \ \ DisplayFunction\  -> \ 
                Identity]; 
functiontable\  = \ \ Table[{\(1\/\(2\ dX\^2\)\) \((\((Y[[2\ k + 
                                1]] - 
                            2\ Y[[2\ k]] + Y[[2\ k - 
                                1]])\)\ \((x - X[[2\ k]])\)\^2)\) + \(\((Y[[2\ k + 1]] - Y[[2\ k - 1]])\)\ \((x - X[[2\ k]])\)\)\/\(2\ dX\) + 
                  Y[[2\ k]], 
                X[[2\ k - 
                        1]] ≤ x && 
                  x ≤ 
                    X[[2\ k + 1]]}, {k, 
                1, n\/2}]; 
g[x_]\  := \ \ Piecewise[
              functiontable]; 
plot2\  = \ 
            FilledPlot[g[x], \ \ {x, \ a, \ b}, 
              Fills -> \ {GrayLevel[ .8]}, \ 
PlotRange -> \ All, \ \ DisplayFunction\  -> \ 
                Identity]; 
plot3\  = \ 
            Graphics[
              verticals]; 

simpsonsum\  = \ N[\((dX/3)\) \((Y[\([1]\)]\  + \ Y[\([n + 1]\)]\  + \ 
                    4  Sum[Y[\([2  k]\)], \ {k, \ 1, \ n/2}]\ \  + \ 
                    2\ Sum[Y[\([2  k + 1]\)], \ {k, \ 1, \ 
                          n/2 - 1}])\)]; 
label\  = \ 
            simpsonsum\ \ "  =  Simpson  "\ \ ; \ 
Show[\ plot2, \ \ plot3, \ 
            plot1, 
PlotLabel\  -> \ 
              label, \ 
\ 
            DisplayFunction\  -> \ $DisplayFunction]
] /; \ EvenQ[n]\)\)\)

Numerical rules

LeftSum[n_] := 
      Module[{X,  Y}, 
X  =  
          Table[a  +  k*(b - a)/n,  {k,  0,  
              n}]; 
Y  =    Map[f,  
            X]; 
 
        N[Sum[Y[[k]],   {k,  1,  n}] *(b - a)/
              n]
]
MidpointSum[n_] := 
      Module[{X,  Y,  midpoints}, 
X  =  
          Table[a  +  k*(b - a)/n,  {k,  0,  
              n}]; 
midpoints  =  
          Table[(X[[k]]  +  X[[k + 1]])/2,   {k,  1,  
              n}]; 
Y  =    Map[f,  
            midpoints]; 
N[
          Sum[Y[[k]],   {k,  1,  n}] *(b - a)/
              n]
]
RightSum[n_] := 
      Module[{X,  Y}, 
X  =  
          Table[a  +  k*(b - a)/n,  {k,  0,  
              n}]; 
Y  =    Map[f,  
            X]; 
 
        N[Sum[Y[[k + 1]],   {k,  1,  n}] *(b - a)/
              n]
]
TrapezoidSum[n_] := 
      Module[{X,  Y}, 
X  =  
          Table[a  +  k*(b - a)/n,  {k,  0,  
              n}]; 
Y  =    Map[f,  
            X]; 
( .5*Y[[1]]  +   .5*
                Y[[n + 1]] +  
              N[Sum[Y[[k]],   {k,  2,  n}] ])*(b - a)/
            n
]
SimpsonSum[
        n_]  :=   Module[{X,  Y,  dX}, 
X  =  
            Table[a  +  k*(b - a)/n,  {k,  0,  
                n}]; 
Y  =    Map[f,  
              X]; 
dX  =  (b - a)/
              n; 
N[(dX/3) (Y[[1]]  +  
                  Y[[n  + 1]]  +  
                  4  Sum[Y[[2  k]],  {k,  1,  n/2}]   +  
                  2 Sum[Y[[2  k + 1]],  {k,  1,  
                        n/2 - 1}])]
] /;  
        EvenQ[n]

f[x_] :=  x^3 Exp[-x^2]
a = 1;
b = 4;
Plot[f[x],  {x,  a,  b}]
functionplot = 
        Plot[f[x], {x, a, b}, DisplayFunction -> Identity,  
          PlotRange ->  {0,  1/2}];

Few intervals

n = 100

    
trapezoid[n]
simpson[n]

\(\(\(
\)\(\(T\  = \ TrapezoidSum[n];\)

    \(S\  = \ SimpsonSum[n];\)

    Print[TableForm[Table[{T, \ S}]]]\)\)\)

Many intervals

SimpsonSum[200]
SimpsonSum[500]

NIntegrate[f[x],  {x,  1,  4}]

Computation for error estimate

k[x_]  =  (∂_x) (∂_x) (∂_x) (∂_x) f[x]

Plot[k[x],  {x,  1,  4}]
Plot[k[x],  {x, 1,  1.2}]

    
FindRoot[k'[x] == 0,  {x,  1.0}]
k[x] /. %

FindRoot[k'[x] == 0,  {x,  2.0}]
k[x] /. %

k[1.07]

estimate[
        n_] := ((20.0 (((b - a))^5))/(180 (n^4)))
Clear[n]
K  =  20;

    
n  = 
      N[ ((((K (((b - a))^5))/(180 )) 10^5)^(1/4))]
estimate[42]

It suffices to take n = 42 to achieve error less than (1/(10^5))  with Simpson's rule.