import numpy as np
# n = number of intervals
# There are n+1 points
def trapz(a,b,n,f,df):
h = (b-a)/n
x = np.linspace(a,b,n+1)
y = f(x)
res = np.sum(y[1:n]) + 0.5*(y[0] + y[n])
return h*res, h*res - (h**2/12)*(df(b) - df(a))
The exact integral is $-\frac{1}{2}(1+\exp(\pi))$.
f = lambda x: np.exp(x)*np.cos(x)
df = lambda x: np.exp(x)*(np.cos(x) - np.sin(x))
qe = -0.5*(1.0 + np.exp(np.pi)) # Exact integral
n,N = 4,10
e1,e2 = np.zeros(N),np.zeros(N)
for i in range(N):
e1[i],e2[i] = trapz(0.0,np.pi,n,f,df) - qe
if i > 0:
print('%6d %24.14e %10.5f %24.14e %10.5f'%(n,e1[i],e1[i-1]/e1[i],e2[i],e2[i-1]/e2[i]))
else:
print('%6d %24.14e %10.5f %24.14e %10.5f'%(n,e1[i],0,e2[i],0))
n = 2*n
4 -1.26567653098185e+00 0.00000 -2.47437900765224e-02 0.00000 8 -3.11816113365945e-01 4.05905 -1.58292813961225e-03 15.63166 16 -7.76577835071954e-02 4.01526 -9.94872006128134e-05 15.91087 32 -1.93958006245669e-02 4.00385 -6.22654792081789e-06 15.97791 64 -4.84778281250620e-03 4.00096 -3.89293344227326e-07 15.99449 128 -1.21187271271594e-03 4.00024 -2.43329250082525e-08 15.99863 256 -3.02963615787633e-04 4.00006 -1.52084034255040e-09 15.99966 512 -7.57406187865683e-05 4.00002 -9.50493017626286e-11 16.00054 1024 -1.89351368735657e-05 4.00000 -5.94013727095444e-12 16.00120 2048 -4.73378310594796e-06 4.00000 -3.73034936274053e-13 15.92381