# Error in Simpson's 3/8 rule is higher than that of Simpson's 1/3 rule

For a given function $$f(x)$$, I have tried to find its numerical integral using Simpson's 1/3 and Simpson's 3/8 rules.

I then compare the solution from the numerical quadratures to the analytical integral to calculate the errors. I find that the error for the 3/8 rule is double that of the 1/3 rule, while the opposite is expected. Where did I go wrong?

The Python3 function for the 1/3 rule:

def simpsons13(a, b, N):
"""
Calculates the numerical integral of a function f(x) using the Simpson's 1/3rd rule:

F(x) = Σ(0 to (N-2)/2) Δx/3 * (f(x(2i)) + 4f(x(2i + 1)) + f(x(2i + 2)))

Parameters:

a:      The lower limit of the definite integral (real)
b:      The upper limit of the definite integral (real)
N:      A positive, even integer to denote the number of intervals of the function for the integral

Returns I, the numerical integral calculated through Simpson's 1/3rd rule
"""

if N%2 != 0:
return "N must be even"

dx = (b - a)/N
narr = np.linspace(a, b, N+1)   # N intervals corresponds to N + 1 points
I = 0

for i in range(int((N - 2) / 2) + 1):

I = I + dx/3*(f(narr[2*i]) + 4*f(narr[2*i + 1]) + f(narr[2*i + 2]))

return I


The Python3 function for the 3/8 rule:

def simpsons38(a, b, N):
"""
Calculates the numerical integral of a function f(x) using the Simpson's 3/8 rule:

F(x) = Σ(0 to (N-3)/3) 3Δx/8 * (f(x(3i)) + 3f(x(3i + 1)) + 3f(x(3i + 2)) + f(x(3i + 3)))

Parameters:

a:      The lower limit of the definite integral (real)
b:      The upper limit of the definite integral (real)
N:      A positive, even integer to denote the number of intervals of the function for the integral

Returns I, the numerical integral calculated through Simpson's 3/8 rule
"""

if N%3 != 0:
return "N must be divisible by 3"

dx = (b - a)/N
narr = np.linspace(a, b, N+1)
I = 0

for i in range(int((N-3)/3) + 1):

I = I + 3*dx/8 * (f(narr[3*i]) + 3*f(narr[3*i + 1]) + 3*f(narr[3*i + 2]) + f(narr[3*i + 3]))

return I



I consider the simple function $$f(x) = xe^x$$ and find its numerical integral from $$[-\pi, \pi]$$ by considering $$N = 24$$ intervals. I calculate the errors between the numerical and analytical solutions as:

error = np.abs(ans - approx)


I get the error by the 1/3 method to be $$0.003669436$$ and by the 3/8 method to be $$0.00816864$$. The 3/8 rule has more than double the error of the 1/3 rule. Why is this happening?

• Could you check if the number of sub-intervals resp. step size in the error formulas refers to the actual sample points or to segments with internal points. That is, in the latter view you would have applied 1/3 to 12 segments and 3/8 to 8 segments. /// Comparing for equal numbers of function evaluations is of course a fair measure. Aug 16 at 8:33

The simple segment-$$[a,b]$$-with-midpoints error formulas are for the 1/3 rule $$E=-\frac{(b-a)^5}{90·2^5}f^{(4)}(\zeta)$$ and for the 3/8 rule $$E=-\frac{(b-a)^5}{90·72}f^{(4)}(\zeta)$$ Now the first rule has 2 sub-intervals while the second has 3. To get comparable values for equal numbers of function evaluations apply this to an interval with 6 sub-intervals, giving 3 segments for the 1/3 rule and 2 segments for the 3/8 rule. Then the composite 1/3 rule has an error of $$E=-3·\frac{((b-a)/3)^5}{90·2^5}f^{(4)}(\zeta) =-\frac{(b-a)^5}{90·2^5·3^4}f^{(4)}(\zeta) =-\frac{((b-a)/6)^5}{30}f^{(4)}(\zeta)$$ and the composite 3/8 rule $$E=-2·\frac{((b-a)/2)^5}{90·72}f^{(4)}(\zeta) =-\frac{(b-a)^5}{40·2^5·3^4} =-\frac{((b-a)/6)^5}{40/3}f^{(4)}(\zeta)$$ The midpoints $$\zeta$$ can and will be different for each formula.
In total the error of the 3/8 rule is about $$\frac94$$ of the error of the 1/3 rule for the same number of sampling points. This is also exactly what you observed, as $$36·\frac94=81$$.