import numpy as np
The next function implements Trapezoid method.
# Performs Trapezoid quadrature
def integrate(a,b,n,f):
h = (b-a)/n
x = np.linspace(a,b,n+1)
y = f(x)
res1 = 0.5*y[0] + np.sum(y[1:n]) + 0.5*y[n]
return h*res1
We will integrate $$ f(x) = x \sqrt{x}, \qquad x \in [0,1] $$ The exact integral is $2/5$ but $f''$ is not finite. Now we perform convergence test. The Peano kernel error formula gives the error estimate $$ |E_n(f)| \le \frac{3 h^2}{16} $$
f = lambda x: x * np.sqrt(x)
a, b = 0.0, 1.0
Ie = 2.0/5.0 # Exact integral
n, N = 2, 10
e = np.zeros(N)
for i in range(N):
h = (b-a)/n
I = integrate(a,b,n,f)
e[i] = Ie - I
if i > 0:
print('%6d %18.8e %14.5g %18.8e'% (n,e[i],e[i-1]/e[i],3*h**2/16))
else:
print('%6d %18.8e %14.5g %18.8e'%(n,e[i],0,3*h**2/16))
n = 2*n
2 -2.67766953e-02 0 4.68750000e-02 4 -7.01811086e-03 3.8154 1.17187500e-02 8 -1.81246480e-03 3.8721 2.92968750e-03 16 -4.63401302e-04 3.9112 7.32421875e-04 32 -1.17671210e-04 3.9381 1.83105469e-04 64 -2.97398625e-05 3.9567 4.57763672e-05 128 -7.49190899e-06 3.9696 1.14440918e-05 256 -1.88304417e-06 3.9786 2.86102295e-06 512 -4.72540682e-07 3.9849 7.15255737e-07 1024 -1.18449772e-07 3.9894 1.78813934e-07
We still get $O(h^2)$ convergence and the error is smaller than the error estimate.