Newton-Cotes rules for integration

• We briefly discuss the simplest integration of a function from an evaluation of the integrand at equally spaced points.

• We want to integrate $f(x)$ in the interval $[a,b]$ with $N$ points, $\Delta = (b-a)/N$.

 

%matplotlib inline
%config InlineBackend.figure_format = 'retina'
import matplotlib.pyplot as plt
import numpy as np
import scipy.special
plt.style.use('style.mpl')

Rectangular

\begin{equation*} \int_a^b f(x) dx = \Delta \sum_{j=0}^{N-1} f \left( a + j \Delta \right) + \mathcal{O}\left(\frac{1}{N}\right). \end{equation*}

Trapezoidal

\begin{equation*} \int_a^b f(x) dx = \Delta \left [\frac{1}{2}f(a) + \sum_{j=1}^{N-1} f\left(a + j \Delta \right) +\frac{1}{2}f(b)\right] + \mathcal{O}\left(\frac{1}{N^2}\right). \end{equation*}

Simpson

\begin{equation*} \int_a^b f(x) dx = \frac{\Delta}{3} \left \{ f(a) + \sum_{j=1}^{N-1}\left[3-(-1)^j\right] f\left(a + j \Delta \right) +f(b)\right\} + \mathcal{O}\left(\frac{1}{N^4}\right). \end{equation*}

Here is an implementation of these schemes

def quadrature(f, a, b, N, *params, method='simpson'):
    '''Newton-Cotes methods for numerical integration of f(x) over x.'''
    
    x, Δx = np.linspace(a, b, N, retstep=True)
    fx = np.array([f(v, *params) for v in x])
    
    if method == 'rectangular':
        I = np.sum(fx)
        
    if method == 'trapezoidal':
        I = 0.5 * (fx[0] + fx[N-1])
        I += np.sum(fx[1:N-1])
        
    elif method == 'simpson':
        assert(N%2 == 1)
        I = (fx[0] + fx[N-1]) / 3.0
        I += (4/3) * np.sum(fx[1:N:2])
        I += (2/3) * np.sum(fx[2:N-1:2])
        
    return Δx * I

Let's use Simpson's rule to evaluate the $\mathrm{erf}(1/2)$ where

\begin{equation*} \mathrm{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x \mathrm{e}^{-t^2} dt \end{equation*}

def erf_kernel(t):
    '''The error function kernel'''
    return  (2.0/np.sqrt(np.pi)) * np.exp(-t**2)

# check quadrature rules
res = quadrature(erf_kernel, 0, 0.5, 2**12+1, method='simpson')
exact = scipy.special.erf(0.5)

print(f'Simpson -> {res:18.16f}')
print(f'Scipy   -> {exact:18.16}')

Simpson -> 0.5204998778130465
Scipy   -> 0.5204998778130465

Error in Newton-Cotes methods

• Difference between exact result and quadrature with $N$ points for \begin{equation*} \frac{2}{\sqrt{\pi}} \int_0^{1/2} \mathrm{e}^{-t^2} dt \end{equation*}

What happens in higher dimensions?

• A method of order $k$ has an error $\sim 1 / N^k$.

• For example, Simpson's method is order $k=4$.

• This is OK in 1 dimension. But what about $d$ dimensions?

  \begin{equation*} I = \int_{\Omega} d^d x\, f(\mathbf{x}) = \int dx_1 \int dx_2 \cdots \int dx_d\, f(\mathbf{x}). \end{equation*}

• We need $M$ points in every dimension for a total of $N = M^d$ points.

• The error $\sim 1/M^k$.

For an order-$k$ scheme in $d$ dimensions, the error $\sim 1 / N^{k/d}$. This is not great...

What about a Monte Carlo approach?

• Let's take our pebbles: every throw is described by a random variable $x_i$

  \begin{equation*} \langle x_i \rangle = \frac{\pi}{4} \cdot 4 + \left( 1 - \frac{\pi}{4} \right) \cdot 0 = \pi \qquad \langle x^2_i \rangle = 4\pi \qquad \sigma = \sqrt{4\pi-\pi^2} \end{equation*}

• Empirical average $X$

  \begin{equation*} X = \frac{1}{N} \sum_{i=1}^N x_i \qquad \langle X \rangle = \frac{1}{N} \sum_{i=1}^N \langle x_i \rangle = \pi \end{equation*}

• What is the standart deviation on the empirical average?

  \begin{equation*} \langle X^2 \rangle - \langle X \rangle^2 = \frac{1}{N^2} \sum_{i,j=1}^N \langle x_i x_j \rangle - \frac{1}{N^2} \sum_{i,j} \langle x_i \rangle \langle x_j \rangle = \frac{1}{N} \sigma^2 \Rightarrow \quad \boxed{\sigma_X = \frac{\sigma}{\sqrt{N}}} \end{equation*}

Monte Carlo approaches have an error $\sim 1 / \sqrt{N}$ in any dimension!

Outline

• Introduction to Monte Carlo

• Newton–Cotes quadrature

• Importance sampling

• Direct sampling methods

• Markov chain sampling and balance condition

• Metropolis-Hastings algorithm

• The two-dimensional Ising model

• Error bar analysis

• References