import numpy as np
def simpson(a,b,n,f,df3):
h = (b-a)/n
x = np.linspace(a,b,n+1)
y = f(x)
res = 4.0*np.sum(y[1:n:2]) + 2.0*np.sum(y[2:n-1:2]) + y[0] + y[n]
est = -(h**4/180.0)*(df3(b) - df3(a))
return (h/3.0)*res, est
The exact integral is $-\frac{1}{2}(1+\exp(\pi))$.
# Function f
f = lambda x: np.exp(x)*np.cos(x)
# Third derivative of f
df3 = lambda x: -2.0*np.exp(x)*(np.cos(x) + np.sin(x))
qe = -0.5*(1.0 + np.exp(np.pi)) # Exact integral
n,N = 2,10
e = np.zeros(N)
for i in range(N):
integral,est = simpson(0.0,np.pi,n,f,df3)
e[i] = qe - integral
if i > 0:
print('%6d %18.8e %14.5f %18.8e'%(n,e[i],e[i-1]/e[i],est))
else:
print('%6d %18.8e %14.5f %18.8e'%(n,e[i],0,est))
n = 2*n
2 -4.77506763e-01 0.00000 -1.63300203e+00 4 -8.54022966e-02 5.59126 -1.02062627e-01 8 -6.13735917e-03 13.91515 -6.37891419e-03 16 -3.94993112e-04 15.53789 -3.98682137e-04 32 -2.48603363e-05 15.88849 -2.49176335e-05 64 -1.55645818e-06 15.97238 -1.55735210e-06 128 -9.73205481e-08 15.99311 -9.73345060e-08 256 -6.08319262e-09 15.99827 -6.08340663e-09 512 -3.80214971e-10 15.99935 -3.80212914e-10 1024 -2.37658782e-11 15.99836 -2.37633071e-11
The convergence rate is $O(h^4)$ and the error estimate becomes very good as $n$ increases.