randomized_integration.py

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import trapezoid, simpson, quad


repo = "https://raw.githubusercontent.com/nicoguaro/matplotlib_styles/master"
style = repo + "/styles/clean.mplstyle"
plt.style.use(style)
plt.rcParams["font.size"] = 20
plt.rcParams["lines.linewidth"] = 3


np.random.seed(69)

a = 1
b = 20
fun = lambda x: np.cos(x)**2
fun = lambda x: np.sin(x)/x

inte_ex, _ = quad(fun, a, b)

npts_list = np.logspace(1, 6, 50, dtype=int)
intes = []
errors = []
intes_r = []
errors_r  = []
intes_t = []
errors_t = []
intes_s = []
errors_s = []
for npts in npts_list:
    x = np.random.uniform(a, b, npts)
    x = np.sort(x)
    xm = (x[:-1] + x[1:])/2
    dx = xm[1:] - xm[:-1]
    f = fun(x)
    
    # Monte Carlo
    inte = (b - a)*np.mean(f)
    intes.append(inte)
    errors.append(np.abs(inte - inte_ex)/inte_ex)
    
    # Riemman
    inte = np.dot(f[1:-1], dx) + f[0] * (xm[0] - a) + f[-1]*(b - xm[-1])
    intes_r.append(inte)
    errors_r.append(np.abs(inte - inte_ex)/inte_ex)
    
    # Trapezoid
    x2 = np.concatenate(([a], x, [b]))
    f2 = fun(x2)
    inte = trapezoid(f2, x2)
    intes_t.append(inte)
    errors_t.append(np.abs(inte - inte_ex)/inte_ex)
    
    # Simpson
    inte = simpson(f2, x2)
    intes_s.append(inte)
    errors_s.append(np.abs(inte - inte_ex)/inte_ex)
    

plt.figure()
plt.loglog(npts_list, errors, "o", label="Monte Carlo")
plt.loglog(npts_list, errors_r, "^", label="Random midpoint")
plt.loglog(npts_list, errors_t, "v", label="Random trapezoid")
plt.loglog(npts_list, errors_s, "s", label="Random Simpson")
#plt.loglog(npts_list, 1/np.sqrt(npts_list), label="$1/\sqrt{n}$")
#plt.loglog(npts_list, 1/npts_list, label="$1/n$")
#plt.loglog(npts_list, 1/npts_list**2, label="$1/n^2$")
#plt.loglog(npts_list, 1/npts_list**3, label="$1/n^3$")
plt.xlabel("Number of points")
plt.ylabel("Relative error")
plt.legend(loc=3)
plt.show()