Simpson's Rule and Variants
Daniel J. Velleman, Amherst College, Amherst, MA 01002

This notebook contains definitions of commands that will compute Simpson's rule and several variants that are discussed in my paper "Simpson Symmetrized and Surpassed".  These definitions are all in initialization cells, which must be evaluated before the commands can be used.  Once the initialization cells have been evaluated, you can use the following commands:

Simpson[f[x],{x,a,b},n]:  This command uses Simpson's rule to estimate (∫_a^b)f(x)dx, using n subintervals.  It first shows a graph in which the curve y = f(x) is represented by a solid line, the parabolic segments that are used to approximate the graph of f are represented by dashed lines, the sample points are marked with dots, and the area that is used to approximate the integral is shaded.  Then it outputs the numerical value of the Simpson's rule approximation.  There are two options that can be used to change the behavior of the command:

Picture→False will suppress the display of the graph, and will just give the numerical value of the approximation.
Digits→(number of digits) specifies the number of digits of precision to be used in all calculations.

For example:

Simpson[Sin[x],{x,0,7Pi},10] will use Simpson's rule with n = 10 to approximate (∫_0^(7  π))sinx dx, first displaying a picture of the approximation.
Simpson[Sin[x],{x,0,7Pi},10,Picture→False] will do the same thing, but without the picture.
Simpson[Sin[x],{x,0,7Pi},10,Digits→20] will do the same thing, using 20 digits of precision in all calculations.
Simpson[Sin[x],{x,0,7Pi},10,Picture→False,Digits→20] will do the same thing, but without the picture, and using 20 digits of precision.

CentSimp[f[x],{x,a,b},n]:  This is like the Simpson command, except that it uses the central Simpson's rule.

ExpCentSimp[f[x],{x,a,b},n,r]:  This is like the Simpson command, except that it uses the r-expanded central Simpson's rule.  The additional parameter r specifies the ratio of the width of the first and last intervals to the width of all other intervals.  A good choice is r = 2/3.

SqueezedSimpson[f[x],{x,a,b},n,r]:  A bonus: this is a rule that is not discussed in my paper.  This command is the same as the Simpson command, except that the first and last intervals can be a different width from the other intervals.  The additional parameter r specifies the ratio of the width of the first and last intervals to the width of all other intervals.  Thus, r = 1 is equivalent to the ordinary Simpson's rule.  You can often get better accuracy using r = .95.

PetersMaley[f[x],{x,a,b},n]:  This is like the Simpson command, except that it uses the rule proposed by G. O. Peters and C. E. Maley in "Numerical integration over any numbers of equal intervals" (Amer. Math. Monthly 75 (1968), 741–744).  One way to think about this rule is as follows:  For 1 ≤ i ≤ n - 2, approximate f on the interval [(Subscript[x, i]), (Subscript[x, i+1])] with the average of the unique quadratic polynomial that agrees with f at (Subscript[x, i-1]), (Subscript[x, i]), and (Subscript[x, i+1]), and the unique quadratic that agrees with f at (Subscript[x, i]), (Subscript[x, i+1]), and (Subscript[x, i+2]).  On the interval [(Subscript[x, 0]), (Subscript[x, 1])], approximate f with the unique quadratic that agrees with f at (Subscript[x, 0]), (Subscript[x, 1]), and (Subscript[x, 2]), and on the interval [(Subscript[x, n-1]), (Subscript[x, n])], approximate f with the unique quadratic that agrees with f at (Subscript[x, n-2]), (Subscript[x, n-1]), and (Subscript[x, n]).  On each interval, approximate the integral of f with the integral of the approximating quadratic polynomial, and sum these approximations to get an approximation of (∫_a^b)f(x) dx.  The picture displayed by the PetersMaley command is based on this way of thinking about the rule.

There are also commands to compute error bounds for some of these rules.  Note that in each of the following commands, M should be an upper bound on the absolute value of the fourth derivative of f on the interval [(Subscript[x, 0]), (Subscript[x, n])].  In the case of the central Simpson's rule, this interval will be larger than [a, b], and the bound on |(f^( (4)))(x)| must be valid on this larger interval.

SimpsonError[a,b,n,M]:  This gives a bound on the magnitude of the error in using Simpson's rule with n subintervals to estimate (∫_a^b)f(x) dx, if |(f^( (4)))(x)| ≤ M.

CentSimpError[a,b,n,M]:  As above, but for the central Simpson's rule.

ExpCentSimp23Error[a,b,n,M]:  As above, but for the 2/3-expanded central Simpson's rule.  Note that this command is only for the 2/3-expanded rule; it does not apply to r ≠ 2/3.

PetersMaleyError[a,b,n,M]:  As above, but for the Peters-Maley rule.

Finally, there is a command for computing the Peano kernel for any of our rules.  According to Peano's theorem, if K is the Peano kernel for one of our rules and (f^( (4))) is continuous on [(Subscript[x, 0]), (Subscript[x, n])], then the error in using the rule to estimate (∫_a^b)f(x)dx is (∫_(Subscript[x, 0])^(Subscript[x, n]))(f^( (4)))(x)K(x)dx.  Note that in the case of the central Simpson's rule, the interval [(Subscript[x, 0]), (Subscript[x, n])] will be larger than [a, b], and you must integrate over this larger interval to get the correct error.

PeanoKernel[rule,a,b,n,(additional parameters)][x]:  This command finds the Peano kernel for the specified rule applied on the interval [a, b], using n subintervals.  Note that the Peano kernel is a function of x.  rule can be any of Simpson, CentSimp, ExpCentSimp, SqueezedSimpson, or PetersMaley.  For ExpCentSimp and SqueezedSimpson you must specify the value of r as an additional parameter after n.  For all of the rules, you can also add, after all other parameters, Digits→(number of digits) to specify how many digits of precision should be used in calculations.

For example:

Plot[Evaluate[PeanoKernel[ExpCentSimp,0,1,10,2/3][x]],{x,0,1}] will draw a graph of the Peano kernel for the 2/3-expanded central Simpson's rule on the interval [0, 1], with n = 10.  The use of Evaluate speeds things up by computing the Peano kernel once at the beginning, rather than recomputing it for each point plotted.
NIntegrate[Derivative[4][Sin][x] PeanoKernel[Simpson,0,7Pi,10][x],{x,0,7Pi}] will numerically integrate the fourth derivative of the sin function times the Peano kernel for Simpson's rule with n = 10 on the interval [0, 7π].  According to Peano's theorem, this should be equal to the error in using Simpson's rule with n = 10 to estimate (∫_0^(7  π))sinx dx.

Definitions of Commands
Open this group of cells if you want to see how the commands are defined.  If you just want to use the commands, you don't need to open this group of cells.  Just be sure to evaluate all initialization cells before using any commands.

IntRuleDefaults = {Picture -> True, 
          Digits -> $MachinePrecision};

PlotPts[f_, start_, end_, numpts_, includeOutline_] := 
      Module[{curve}, 
        curve = Cases[
              Flatten[Plot[f[x], {x, start, end}, 
                    PlotPoints -> Max[3, numpts], 
                    DisplayFunction -> 
                      Identity][[1]]], _Line][[1, 1]]; 
        If[includeOutline, 
          Join[Drop[
              curve, -1], {{end, f[end]}, {end, 0}, {start, 0}, {start, 
                f[start]}}], curve]]

ParabPts[{x0_, y0_}, {x1_, y1_}, {x2_, y2_}, start_, end_, numpts_, 
        includeOutline_] := 
      Module[{A, B, C, parab}, 
        parab[x_] := 
          A (x^2) + B x + C; {A, B, 
            C} = {A, B, 
              C} /. Solve[{parab[x0] == y0, 
                  parab[x1] == y1, parab[x2] == y2}, {A, B, 
                  C}][[1]]; 
        PlotPts[parab, start, end, numpts, includeOutline]]

ParabPts[y0_, y1_, y2_, x0_, h_, start_, end_, numpts_, 
        includeOutline_] := 
      Module[{A, B, parab}, A = ((y0 - 2 y1 + y2)/(2 (h^2))); 
        B = ((-3 y0 + 4 y1 - y2)/(2 h)); 
        parab[x_] := A (((x - x0))^2) + B (x - x0) + y0; 
        PlotPts[parab, start, end, numpts, includeOutline]]

ParabPts[y0_, y1_, y2_, y3_, x0_, h_, start_, end_, numpts_, 
        includeOutline_] := 
      Module[{A, B, C, parab}, A = ((y0 - y1 - y2 + y3)/(4 (h^2))); 
        B = ((-3 y0 - y1 + 7  y2 - 3  y3)/(4 h)); 
        C = ((y0 + 3 y1 - 3 y2 + y3)/2); 
        parab[x_] := A (((x - x0))^2) + B (x - x0) + C; 
        PlotPts[parab, start, end, numpts, includeOutline]]

Outline[pts_] := {AbsoluteDashing[{4}], Line[Drop[pts, -3]], 
        AbsoluteDashing[{}], Line[Take[pts, -4]]}

\(Simpson[funcexp_, {var_, a_, b_}, n_Integer, opts___] := 
      Module[{f, p, dx, xx, yy, pts, pic, 
          dig}, {pic, 
            dig} = \({Picture, 
                Digits} /. {opts}\) /. IntRuleDefaults; 
        If[EvenQ[n] && n ≥ 2, f = Function[var, funcexp]; 
          dx = N[\(b - a\)\/n, dig]; 
          xx = Table[N[a + dx\ i, dig], {i, 0, n}]; 
          yy = Table[
              N[f[xx[[i]]], dig], {i, 1, 
                n + 1}]; 
          If[pic, p = 
              Plot[funcexp, {var, a, b}, DisplayFunction -> Identity]; 
            pts = Join[
                Join @@ Table[
                    Drop[ParabPts[
                        yy[[i]], 
                        yy[[i + 1]], 
                        yy[[i + 2]], 
                        xx[[i]], dx, 
                        xx[[i]], 
                        xx[[i + 2]], 
                        Ceiling[50\/n], False], \(-1\)], {i, 1, n - 1, 
                      2}], {{b, 
                    yy[[n + 1]]}, {b, 
                    0}, {a, 0}, {a, 
                    yy[[1]]}}]; 
            Show[{Graphics[{GrayLevel[ .8], Polygon[pts], GrayLevel[0], 
                    Outline[pts], 
                    Table[Line[{{xx[[
                              i]], 
                            0}, {xx[[
                              i]], 
                            yy[[
                              i]]}}], {i, 3, n - 1, 2}], 
                    AbsolutePointSize[5], 
                    Table[Point[{xx[[
                            i]], 
                          yy[[i]]}], {i, 
                        1, n + 1}]}], p}, Axes -> True, 
              PlotRange -> All]]; 
          1\/3\ dx\ \((yy[[1]] + 
                4\ yy[[n]] + 
                yy[[n + 1]] + 
                2\ Sum[2\ yy[[i]] + 
                      yy[[i + 1]], {i, 
                      2, n - 2, 2}])\), 
          Print["n must be even and at least 2 for Simpson's rule"]]]; 
    SetAttributes[Simpson, HoldAll]\)

\(CentSimp[funcexp_, {var_, a_, b_}, n_Integer, opts___] := 
      Module[{f, p, dx, xx, yy, pts, pic, 
          dig}, {pic, 
            dig} = \({Picture, 
                Digits} /. {opts}\) /. IntRuleDefaults; 
        If[n ≥ 3, f = Function[var, funcexp]; 
          dx = N[\(b - a\)\/\(n - 1\), dig]; 
          xx = Table[N[a - dx\/2 + dx\ i, dig], {i, 0, n}]; 
          yy = Table[
              N[f[xx[[i]]], dig], {i, 1, 
                n + 1}]; 
          If[pic, p = 
              Plot[funcexp, {var, 
                  xx[[1]], 
                  xx[[n + 1]]}, 
                DisplayFunction -> Identity]; 
            pts = Table[
                ParabPts[yy[[i - 1]], 
                  yy[[i]], 
                  yy[[i + 1]], 
                  xx[[i - 1]], dx, 
                  xx[[i]] - dx\/2, 
                  xx[[i]] + dx\/2, 
                  Ceiling[25\/n], True], {i, 2, n}]; 
            Show[{Graphics[{GrayLevel[ .8], 
                    Table[Polygon[
                        pts[[i]]], {i, 
                        1, n - 1}], GrayLevel[0], 
                    Table[Outline[
                        pts[[i]]], {i, 
                        1, n - 1}], AbsolutePointSize[5], 
                    Table[Point[{xx[[
                            i]], 
                          yy[[i]]}], {i, 
                        1, n + 1}]}], p}, Axes -> True, 
              PlotRange -> All]]; 
          dx\ \((1\/24\ \((yy[[1]] + 
                      yy[[n + 1]] + 
                      23\ \((yy[[2]] + 
                            yy[[
                              n]])\))\) + ∑\+\(i = 3\)\%\(n - 1\)yy[[i]])\), 
          Print["n must be at least 3 for CentSimp."]]]; 
    SetAttributes[CentSimp, HoldAll]\)

\(ExpCentSimp[funcexp_, {var_, a_, b_}, n_Integer, r_, opts___] := 
      Module[{f, p, dx, xx, yy, pts, pic, 
          dig}, {pic, 
            dig} = \({Picture, 
                Digits} /. {opts}\) /. IntRuleDefaults; 
        If[n ≥ 5, f = Function[var, funcexp]; 
          dx = N[\(b - a\)\/\(n - 2 + 2\ r\), dig]; 
          xx = Join[{a}, 
              Table[N[a + \((r + i)\)\ dx, dig], {i, 0, n - 2}], {b}]; 
          yy = Table[
              N[f[xx[[i]]], dig], {i, 1, 
                n + 1}]; 
          If[pic, p = 
              Plot[funcexp, {var, a, b}, DisplayFunction -> Identity]; 
            pts = Join[{ParabPts[{xx[[1]], yy[[1]]}, {xx[[2]], 
                      yy[[2]]}, {xx[[3]], 
                      yy[[3]]}, 
                    xx[[1]], 
                    xx[[2]] + dx\/2, 
                    Ceiling[25\/n], True]}, 
                Table[ParabPts[
                    yy[[i - 1]], 
                    yy[[i]], 
                    yy[[i + 1]], 
                    xx[[i - 1]], dx, 
                    xx[[i]] - dx\/2, 
                    xx[[i]] + dx\/2, 
                    Ceiling[25\/n], True], {i, 3, 
                    n - 1}], {ParabPts[{xx[[
                        n - 1]], 
                      yy[[n - 
                          1]]}, {xx[[
                        n]], 
                      yy[[
                        n]]}, {xx[[
                        n + 1]], 
                      yy[[n + 1]]}, 
                    xx[[n]] - dx\/2, 
                    xx[[n + 1]], 
                    Ceiling[25\/n], True]}]; 
            Show[{Graphics[{GrayLevel[ .8], 
                    Table[Polygon[
                        pts[[i]]], {i, 
                        1, n - 1}], GrayLevel[0], 
                    Table[
                      Outline[
                        pts[[i]]], {i, 
                        1, n - 1}], AbsolutePointSize[5], 
                    Table[Point[{xx[[
                            i]], 
                          yy[[i]]}], {i, 
                        1, n + 1}]}], p}, Axes -> True, 
              PlotRange -> All]]; 
          dx\ \((\(\((2\ r + 1)\)\ \((2\ r\^2 + 2\ r - 1)\)\ \((yy[[1]] + yy[[n + 1]])\)\)\/\(12\ r\ \((r + 1)\)\) + \((\(\((2\ r + 1)\)\^2\ \((r + 2)\)\)\/\(24\ r\) + 
                      1\/24)\)\ \((yy[[2]] + 
                      yy[[n]])\) + \((\(\((2\ r + 1)\)\^2\ \((1 - r)\)\)\/\(24\ \((r + 1)\)\) + 
                      23\/24)\)\ \((yy[[3]] + 
                      yy[[n - 
                          1]])\) + ∑\+\(i = 4\)\%\(n - 2\)yy[[i]])\), 
          Print["n must be at least 5 for ExpCentSimp."]]]; 
    SetAttributes[ExpCentSimp, HoldAll]\)

\(SqueezedSimpson[funcexp_, {var_, a_, b_}, n_Integer, r_, opts___] := 
      Module[{f, p, dx, xx, yy, pts, pic, 
          dig}, {pic, 
            dig} = \({Picture, 
                Digits} /. {opts}\) /. IntRuleDefaults; 
        If[EvenQ[n] && \((n ≥ 6)\), f = Function[var, funcexp]; 
          dx = N[\(b - a\)\/\(n - 2 + 2\ r\), dig]; 
          xx = Join[{a}, 
              Table[N[a + \((r + i)\)\ dx, dig], {i, 0, n - 2}], {b}]; 
          yy = Table[
              N[f[xx[[i]]], dig], {i, 1, 
                n + 1}]; 
          If[pic, p = 
              Plot[funcexp, {var, a, b}, DisplayFunction -> Identity]; 
            pts = Join[
                Drop[ParabPts[{xx[[1]], 
                      yy[[1]]}, {xx[[2]], 
                      yy[[2]]}, {xx[[3]], 
                      yy[[3]]}, 
                    xx[[1]], 
                    xx[[3]], 
                    Ceiling[50\/n], False], \(-1\)], 
                Join @@ Table[
                    Drop[ParabPts[
                        yy[[i]], 
                        yy[[i + 1]], 
                        yy[[i + 2]], 
                        xx[[i]], dx, 
                        xx[[i]], 
                        xx[[i + 2]], 
                        Ceiling[50\/n], False], \(-1\)], {i, 3, n - 3, 2}], 
                Drop[ParabPts[{xx[[
                        n - 1]], 
                      yy[[n - 
                          1]]}, {xx[[
                        n]], 
                      yy[[
                        n]]}, {xx[[
                        n + 1]], 
                      yy[[n + 1]]}, 
                    xx[[n - 1]], 
                    xx[[n + 1]], 
                    Ceiling[50\/n], False], \(-1\)], {{b, 
                    yy[[n + 1]]}, {b, 
                    0}, {a, 0}, {a, 
                    yy[[1]]}}]; 
            Show[{Graphics[{GrayLevel[ .8], Polygon[pts], GrayLevel[0], 
                    Outline[pts], 
                    Table[Line[{{xx[[i]], 0}, {xx[[i]], 
                            yy[[i]]}}], {i, 3, n - 1, 2}], AbsolutePointSize[5], 
                    Table[Point[{xx[[i]], yy[[i]]}], {i, 1, n + 1}]}], p}, 
              Axes -> True, PlotRange -> All]]; 
          dx\ \((\(\((2\ r - 1)\)\ \((r + 1)\)\ \((yy[[1]] + yy[[n + 1]])\)\)\/\(6\ r\) + \(\((r + 1)\)\^3\ \((yy[[2]] + yy[[n]])\)\)\/\(6\ r\) + \((1\/6\ \((2 - r)\)\ \((r + 1)\) + 
                      1\/3)\)\ \((yy[[3]] + 
                      yy[[n - 1]])\) + 
                2\/3\ Sum[
                    2\ yy[[i]] + 
                      yy[[i + 1]], {i, 
                      4, n - 4, 2}] + 
                4\/3\ yy[[n - 2]])\), 
          Print["n must be even and at least 6 for Squeezed Simpson's rule"]]]; SetAttributes[SqueezedSimpson, HoldAll]\)

\(PetersMaley[funcexp_, {var_, a_, b_}, n_Integer, opts___] := 
      Module[{f, p, dx, xx, yy, pts, pic, 
          dig}, {pic, 
            dig} = \({Picture, 
                Digits} /. {opts}\) /. IntRuleDefaults; 
        If[n ≥ 5, f = Function[var, funcexp]; 
          dx = N[\(b - a\)\/n, dig]; 
          xx = Table[N[a + i\ dx, dig], {i, 0, n}]; 
          yy = Table[
              N[f[xx[[i]]], dig], {i, 1, 
                n + 1}]; 
          If[pic, p = 
              Plot[funcexp, {var, a, b}, DisplayFunction -> Identity]; 
            pts = Join[
                Drop[ParabPts[yy[[1]], 
                    yy[[2]], 
                    yy[[3]], 
                    xx[[1]], dx, 
                    xx[[1]], 
                    xx[[2]], 
                    Ceiling[25\/n], False], \(-1\)], 
                Join @@ Table[
                    Drop[ParabPts[
                        yy[[i - 1]], 
                        yy[[i]], 
                        yy[[i + 1]], 
                        yy[[i + 2]], 
                        xx[[i - 1]], dx, 
                        xx[[i]], 
                        xx[[i + 1]], 
                        Ceiling[25\/n], False], \(-1\)], {i, 2, n - 1}], 
                Drop[ParabPts[
                    yy[[n - 1]], 
                    yy[[n]], 
                    yy[[n + 1]], 
                    xx[[n - 1]], dx, 
                    xx[[n]], 
                    xx[[n + 1]], 
                    Ceiling[25\/n], False], \(-1\)], {{b, 
                    yy[[n + 1]]}, {b, 
                    0}, {a, 0}, {a, 
                    yy[[1]]}}]; 
            Show[{Graphics[{GrayLevel[ .8], Polygon[pts], GrayLevel[0], 
                    Outline[pts], 
                    Table[Line[{{xx[[
                              i]], 
                            0}, {xx[[
                              i]], 
                            yy[[
                              i]]}}], {i, 2, n}], 
                    AbsolutePointSize[5], 
                    Table[Point[{xx[[
                            i]], 
                          yy[[i]]}], {i, 
                        1, n + 1}]}], p}, Axes -> True, 
              PlotRange -> All]]; 
          dx\ \((3\/8\ \((yy[[1]] + 
                      yy[[n + 1]])\) + 
                7\/6\ \((yy[[2]] + 
                      yy[[n]])\) + 
                23\/24\ \((yy[[3]] + 
                      yy[[n - 
                          1]])\) + ∑\+\(i = 4\)\%\(n - 2\)yy[[i]])\), 
          Print["n must be at least 5 for Peters-Maley"]]]; 
    SetAttributes[PetersMaley, HoldAll]\)

SimpsonError[a_, b_, n_, 
        M_] := (((((b - a))^5) M)/(180 (n^4)))

CentSimpError[a_, b_, n_, 
        M_] := ((17 (((b - a))^5) M)/(5760 (((n - 1))^4)))

ExpCentSimp23Error[a_, b_, n_, 
        M_] := (((((b - a))^5) M)/(180 (((n - (2/3)))^5))) (0.1494112327604335088340350030749141271008 n - 
            0.03093888526309228169022850227100283534507)

PetersMaleyError[a_, b_, n_, 
        M_] := (((19 n - 30) (((b - a))^5) M)/(720 (n^5)))

PeanoKernel[rule_, a_, b_, n_Integer, opts___][t_] := 
      Module[{x}, 
        (1/6) (If[t < b, 
                (1/4) ((((b - t))^4) - If[t < a, (((a - t))^4), 0]), 
                0] - rule[If[t < x, (((x - t))^3), 0], {x, a, b}, n, opts, 
                Picture -> False])]

Examples

Integrate[Sin[x], {x, 0, 7  Pi}]

Simpson[Sin[x], {x, 0, 7  Pi}, 10]

Simpson[Sin[x], {x, 0, 7  Pi}, 20,  Digits -> 20]

Simpson[Sin[x], {x, 0, 7  Pi}, 50, Picture -> False, 
      Digits -> 20]

CentSimp[Sin[x], {x, 0, 7  Pi}, 10]

CentSimp[Sin[x], {x, 0, 7  Pi}, 50, Picture -> False, 
      Digits -> 20]

ExpCentSimp[Sin[x], {x, 0, 7  Pi}, 10, 2/3]

ExpCentSimp[Sin[x], {x, 0, 7  Pi}, 50, 2/3, Picture -> False, 
      Digits -> 20]

SqueezedSimpson[Sin[x], {x, 0, 7  Pi}, 10,  .95]

SqueezedSimpson[Sin[x], {x, 0, 7  Pi}, 50, 95/100, 
      Picture -> False, Digits -> 20]

PetersMaley[Sin[x], {x, 0, 7  Pi}, 10]

PetersMaley[Sin[x], {x, 0, 7  Pi}, 50, Picture -> False, 
      Digits -> 20]

N[SimpsonError[0, 7  Pi, 20, 1]]

Integrate[Sin[x], {x, 0, 7  Pi}] - 
      Simpson[Sin[x], {x, 0, 7  Pi}, 20, Picture -> False]

N[ExpCentSimp23Error[0, 7  Pi, 20, 1]]

Integrate[Sin[x], {x, 0, 7  Pi}] - 
      ExpCentSimp[Sin[x], {x, 0, 7  Pi}, 20, 2/3, 
        Picture -> False]

SimpsonPK = 
      Plot[Evaluate[PeanoKernel[Simpson, 0, 1, 10][x]], {x, 0, 
          1}]

ExpCentSimpPK = 
      Plot[Evaluate[PeanoKernel[ExpCentSimp, 0, 1, 10, 2/3][x]], {x, 0, 
          1}]

SqueezedSimpsonPK = 
      Plot[Evaluate[PeanoKernel[SqueezedSimpson, 0, 1, 10,  .95][x]], {x, 
          0, 1}]

PetersMaleyPK = 
      Plot[Evaluate[PeanoKernel[PetersMaley, 0, 1, 10][x]], {x, 0, 
          1}]

Show[SimpsonPK, ExpCentSimpPK, SqueezedSimpsonPK, 
      PetersMaleyPK]

NIntegrate[Derivative[4][Sin][
          x] PeanoKernel[Simpson, 0, 7  Pi, 10][x], {x, 0, 
        7  Pi}]

Integrate[Sin[x], {x, 0, 7  Pi}] - 
      Simpson[Sin[x], {x, 0, 7  Pi}, 10, Picture -> False]

NIntegrate[Derivative[4][Sin][
          x] PeanoKernel[ExpCentSimp, 0, 7  Pi, 10, 2/3][x], {x, 0, 
        7  Pi}]

Integrate[Sin[x], {x, 0, 7  Pi}] - 
      ExpCentSimp[Sin[x], {x, 0, 7  Pi}, 10, 2/3, 
        Picture -> False]

NIntegrate[Derivative[4][Sin][
          x] PeanoKernel[CentSimp, 0, 7  Pi, 10][x], {x, -7 Pi/18, 
        7  Pi (19/18)}]

Integrate[Sin[x], {x, 0, 7  Pi}] - 
      CentSimp[Sin[x], {x, 0, 7  Pi}, 10, Picture -> False]