Error Simpson's rule for polynomial not equal to 0

Question

I implemented Simpson's rule in Python and tested the error for a cubic polynomial. The error bound of Simpson's rule states that it should be equal to 0, but it is not. Does anybody know why?

def simps(f,a,b,N=50):
   if N % 2 == 1:
      raise ValueError("N must be an even integer.")
   dx = (b-a)/N
   x = np.linspace(a,b,N+1)
   y = f(x)
   S = dx/3 * np.sum(y[0:-1:2] + 4*y[1::2] + y[2::2])
   return S

I expected the error to be equal to 0 for each $n$, but it was only equal to 0 for small $n$.

Answer 1

Firstly, the error will typically not be 0: as simpson's rule is a numerical approximation, there will usually be an error. Sometimes we just happen to obtain the exact answer, e.g.

f = lambda x: x
F = lambda x: x ** 2 / 2
# Error for default N
simps(f, 0, 1) - (F(1) - F(0))
# 0.0
# Error for custom N
simps(f, 0, 1, 100) - (F(1) - F(0))
# 0.0

Here is an example with exp:

f = np.exp
F = np.exp
# Error for default N
simps(f, 0, 1) - (F(1) - F(0))
# 1.527289184011238e-09
# Error for custom N
simps(f, 0, 1, 100) - (F(1) - F(0))
# 9.545919610332021e-11

As you can see, we obtain quite small error terms, but still are not able to obtain the exact value.

The key takeaway here is twofold:

• Numerical approximations are expected to have errors
• Computer arithmetic introduces a whole new world of errors: a numerical approximation might have on paper an error of magnitude 1e-10 but due to floating point precision that could jump to 1e-6.