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

Question

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?

Answer 1

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$.

---

Note that if this gets expanded to Runge-Kutta methods, where the "classical" method is based on the 1/3 rule and the 3/8 method one the 3/8 rule, both methods have 4 stages, that is, 4 function evaluations. This means that both methods can be compared based on the one-segment error formulas, where thus the 3/8 method has the smaller step error. (One would have to confirm that with the actual one-step Runge-Kutta error, which was indeed done already by Kutta in the original paper.)