import numpy as np
The next function implements Trapezoid and Simpson methods.
# Performs Trapezoid and Simpson 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]
res2 = 4.0*np.sum(y[1:n:2]) + 2.0*np.sum(y[2:n-1:2]) + y[0] + y[n]
return h*res1, (h/3.0)*res2
The next function performs convergence test.
def test(a,b,f,Ie,n,N):
e1,e2 = np.zeros(N), np.zeros(N)
for i in range(N):
I1,I2 = integrate(a,b,n,f)
e1[i],e2[i] = Ie - I1, Ie - I2
if i > 0:
print('%6d %18.8e %14.5g %18.8e %14.5g'%
(n,e1[i],e1[i-1]/e1[i],e2[i],e2[i-1]/e2[i]))
else:
print('%6d %18.8e %14.5g %18.8e %14.5g'%(n,e1[i],0,e2[i],0))
n = 2*n
The exact integral is $2/9$.
f = lambda x: x**3 * np.sqrt(x)
a, b = 0.0, 1.0
Ie = 2.0/9.0 # Exact integral
n, N = 2, 10
test(a,b,f,Ie,n,N)
2 -7.19719516e-02 0 -3.37000954e-03 0 4 -1.81666065e-02 3.9618 -2.31491460e-04 14.558 8 -4.55322411e-03 3.9898 -1.54299774e-05 15.003 16 -1.13906170e-03 3.9973 -1.00755946e-06 15.314 32 -2.84814093e-04 3.9993 -6.48919993e-08 15.527 64 -7.12066288e-05 3.9998 -4.14075488e-09 15.672 128 -1.78018541e-05 4 -2.62556588e-10 15.771 256 -4.45047596e-06 4 -1.65759628e-11 15.84 512 -1.11261977e-06 4 -1.04333209e-12 15.888 1024 -2.78154993e-07 4 -6.55309140e-14 15.921
The exact integral is $\arctan(5-\pi) + \arctan(\pi)$.
f = lambda x: 1.0/(1.0 + (x-np.pi)**2)
a, b = 0.0, 5.0
Ie = np.arctan(b-np.pi) + np.arctan(np.pi)
n, N = 2, 10
test(a,b,f,Ie,n,N)
2 1.73111440e-01 0 -2.85329178e-01 0 4 7.10986397e-02 2.4348 3.70943729e-02 -7.692 8 7.49581780e-03 9.4851 -1.37051228e-02 -2.7066 16 1.95340266e-03 3.8373 1.05930941e-04 -129.38 32 4.89160655e-04 3.9934 1.07998844e-06 98.085 64 1.22340738e-04 3.9983 6.74323792e-08 16.016 128 3.05883472e-05 3.9996 4.21691793e-09 15.991 256 7.64728450e-06 3.9999 2.63594035e-10 15.998 512 1.91183348e-06 4 1.64748215e-11 16 1024 4.77959142e-07 4 1.02939879e-12 16.004
The exact integral is $2/3$.
f = lambda x: np.sqrt(x)
a, b = 0.0, 1.0
Ie = 2.0/3.0
n, N = 2, 10
test(a,b,f,Ie,n,N)
2 6.31132761e-02 0 2.85954792e-02 0 4 2.33836204e-02 2.699 1.01404019e-02 2.82 8 8.53644504e-03 2.7393 3.58738658e-03 2.8267 16 3.08546979e-03 2.7667 1.26847804e-03 2.8281 32 1.10773039e-03 2.7854 4.48483920e-04 2.8284 64 3.95855288e-04 2.7983 1.58563588e-04 2.8284 128 1.41009370e-04 2.8073 5.60607304e-05 2.8284 256 5.01176901e-05 2.8136 1.98204636e-05 2.8284 512 1.77851167e-05 2.818 7.00759224e-06 2.8284 1024 6.30444768e-06 2.821 2.47755801e-06 2.8284
Since $f'(0)$ is not finite, we do not get the optimal convergence rates. But the errors still decrease and both methods converge at same rate.
The exact integral is 7.95492652101284.
f = lambda x: np.exp(np.cos(x))
a, b = 0.0, 2*np.pi
Ie = 7.95492652101284
n, N = 2, 5
test(a,b,f,Ie,n,N)
2 -1.74053505e+00 0 7.20800573e-01 0 4 -3.43969188e-02 50.601 5.34315792e-01 1.349 8 -1.25168894e-06 27480 1.14639707e-02 46.608 16 -3.55271368e-15 3.5232e+08 4.17229640e-07 27476 32 -3.55271368e-15 1 -1.77635684e-15 -2.3488e+08
The integrand is periodic and the error formula suggests that we should expect fast convergence.