Numerical Integration ¶

Contents¶

Define and display integral
sympy.integrate
Midpoint rule
Trapezoidal rule
Gauss–Legendre quadrature
Monte Carlo integration
scipy.integrate.quad
Summary

import integration_helpers as ih
import sympy as sp

Define and display integral ¶

Define $f$ as the function we want to integrate from $a$ to $b$ ($f$ has to be analytically integrable).

x = sp.Symbol('x')
f = sp.sin(x) * x**2
a = 7
b = 15

$$F = \int_{7}^{15} x^{2} \sin{\left(x \right)}\,\mathrm{d}x$$

Plot function and fill out the integral.

fn = sp.lambdify(x, f, modules='numpy')
ih.plot_int(fn, a, b)

Analytical integration via sympy.integrate ¶

$$\int_a^b f(x)\,\mathrm{d}x = F(b) - F(a)$$

Generate a reference value for error calculations.

def int_sympy(f, a, b, x):
    F = integrate(f, x)
    F = lambdify(x, F, modules='numpy')
    return F(b) - F(a)

F_sympy = ih.int_sympy(f, a, b, x)
t_sympy = %timeit -q -o ih.int_sympy(f, a, b, x)
t_sympy = t_sympy.average * 1000

$$F_{\mathrm{sympy}} = 215.154633$$

$$t_{\mathrm{sympy}} = 192.151558\,\mathrm{ms}$$

Numerical integration via midpoint rule ¶

$$\int_a^b f(x)\,\mathrm{d}x \approx \frac{b - a}{n} \sum_{i=1}^n f(m_i)$$

with $m_i$ as the midpoint of the $i$-th interval.

Plot function and rectangles one obtains while using the chained midpoint rule.

n_midpoint = 10
ih.plot_midpoint(fn, a, b, n_midpoint)

def int_midpoint(f, a, b, n):
    width = (b - a) / n
    mi = np.arange(a + width / 2, b, width)
    return width * np.sum(f(mi))

F_midpoint = ih.int_midpoint(fn, a, b, n_midpoint)
e_midpoint = 100 * (F_midpoint - F_sympy) / F_sympy
t_midpoint = %timeit -q -o ih.int_midpoint(fn, a, b, n_midpoint)
t_midpoint = t_midpoint.average * 1000

$$F_{\mathrm{midpoint}} = 220.507528$$

$$\epsilon_{\mathrm{midpoint}} = 2.487929e+00\,\%$$ $$t_{\mathrm{midpoint}} = 0.022549\,\mathrm{ms}$$

Numerical integration via trapezoidal rule ¶

$$\int_a^b f(x)\,\mathrm{d}x \approx \frac{b - a}{2 \cdot n} \sum_{i=1}^n f(x_{i-1}) + f(x_i)$$

Plot function and trapezoids one obtains while using the chained trapezoidal rule.

n_trapez = 10
ih.plot_trapez(fn, a, b, n_trapez)

def int_trapez(f, a, b, n):
    width = (b - a) / n
    xi = np.arange(a + width, b, width)
    return 0.5 * width * (f(a) + np.sum(2 * f(xi)) + f(b))

F_trapez = ih.int_trapez(fn, a, b, n_trapez)
e_trapez = 100 * (F_trapez - F_sympy) / F_sympy
t_trapez = %timeit -q -o ih.int_trapez(fn, a, b, n_trapez)
t_trapez = t_trapez.average * 1000

$$F_{\mathrm{trapez}} = 204.521496$$

$$\epsilon_{\mathrm{trapez}} = -4.942091e+00\,\%$$ $$t_{\mathrm{trapez}} = 0.031264\,\mathrm{ms}$$

Numerical integration via Gauss–Legendre quadrature ¶

$$\int_a^b f(x)\,\mathrm{d}x \approx \frac{b - a}{2} \sum_{i=1}^n \omega_i \cdot f\left(\frac{b - a}{2} \xi_i + \frac{a + b}{2}\right)$$

with $\xi_i$ as the roots of the Legendre-polynomial $$P_n(x) = \sum_{i=0}^{\left \lfloor \frac{n}{2} \right \rfloor} (-1)^k \cdot \frac{(2n-2k)!}{2^n \cdot k! \cdot (n-k)! \cdot (n-2k)!} \cdot x^{n-2k}$$ and the weights $$\omega_i = \frac{2}{(1-\xi_i^2) \cdot P_n^\prime(\xi_i)^2}$$

def int_gauss(f, a, b, n):
    prefactors = [0] * (n + 1)
    for k in range(n // 2 + 1):
        prefactors[2 * k] = (-1)**k * factorial(2 * (n - k)) / \
                            (factorial(n - k) * factorial(n - 2 * k) * factorial(k) * 2**n)
    pol = np.poly1d(prefactors)
    dpol = pol.deriv()
    roots = pol.r
    weights = 2 / ((1 - roots**2) * dpol(roots)**2)
    return 0.5 * (b - a) * np.sum(weights * f(0.5 * (b - a) * roots + 0.5 * (a + b)))

n_gauss = 10
F_gauss = ih.int_gauss(fn, a, b, n_gauss)
e_gauss = 100 * (F_gauss - F_sympy) / F_sympy
t_gauss = %timeit -q -o ih.int_gauss(fn, a, b, n_gauss)
t_gauss = t_gauss.average * 1000

$$F_{\mathrm{gauss}} = 215.154633$$

$$\epsilon_{\mathrm{gauss}} = -4.921084e-10\,\%$$ $$t_{\mathrm{gauss}} = 0.391606\,\mathrm{ms}$$

Numerical integration via Monte Carlo integration ¶

$$\int_a^b f(x)\,\mathrm{d}x \approx \frac{b - a}{n} \sum_{i=1}^n f(x_i)$$

with $x_i$ as randomly sampled values in the interval $[a, b]$.

Plot function and randomly sampled lines (just for demonstration, they will not be the same in int_mc).

n_mc = 1000
ih.plot_mc(fn, a, b, n_mc)

def int_mc(f, a, b, n):
    xi = uniform(a, b, n)
    return np.sum(f(xi)) * (b - a) / n

F_mc = ih.int_mc(fn, a, b, n_mc)
e_mc = 100 * (F_mc - F_sympy) / F_sympy
t_mc = %timeit -q -o ih.int_mc(fn, a, b, n_mc)
t_mc = t_mc.average * 1000

$$F_{\mathrm{MC}} = 226.283774$$

$$\epsilon_{\mathrm{MC}} = 5.172624e+00\,\%$$ $$t_{\mathrm{MC}} = 0.083587\,\mathrm{ms}$$

Numerical integration via scipy.integrate.quad ¶

The documentation states that an integration technique from the Fortran library QUADPACK will be used.

def int_quad(f, a, b, **kwargs):
    F, _ = quad(f, a, b)
    return F

F_quad = ih.int_quad(fn, a, b)
e_quad = 100 * (F_quad - F_sympy) / F_sympy
t_quad = %timeit -q -o ih.int_quad(fn, a, b)
t_quad = t_quad.average * 1000

$$F_{\mathrm{quad}} = 215.154633$$

$$\epsilon_{\mathrm{quad}} = 7.925939e-14\,\%$$ $$t_{\mathrm{quad}} = 0.059027\,\mathrm{ms}$$

Summary ¶

Method	Integral value	Error [%]	Time [ms]
sympy.integrate	215.154633	-	192.151558
Midpoint rule (n=10)	220.507528	2.487929e+00	0.022549
Trapezoidal rule (n=10)	204.521496	-4.942091e+00	0.031264
Gauss-Legendre rule (n=10)	215.154633	-4.921084e-10	0.391606
Monte Carlo integration (n=1000)	226.283774	5.172624e+00	0.083587
scipy.integrate.quad	215.154633	7.925939e-14	0.059027