phys512/integration.md

Integration

Newton-Cotes methods

A related problem to interpolation is when we need to calculate the integral of a function

I = \int_{x_1}^{x_N} f(x) dx

given values of the function at a discrete set of points x_1\dots x_N, which we'll write as f_i\equiv f(x_i). For simplicity here we'll assume that the spacing between points on the grid is constant x_{i+1}-x_i = \Delta x, but it's straighforward to generalize to non-uniform sampling if you need to.

To estimate the value of the integral, we need to model how the function behaves in each interval x_i<x<x_{i+1}. As with interpolation, we can take f(x) in each interval to be a polynomial, including additional terms in the polynomial to increase the accuracy of our approximation.

Rectangle rule

The first approximation is to assume that f(x) is a constant f(x)=f_i in the interval x_i\leq x\leq x_{i+1}. The total integral is then

I\approx \sum_{i=1}^{N-1} f_i \Delta x.

Trapezoidal rule

If we make a linear approximation, we can write

f(x) \approx f_i + {f_{i+1}-f_i\over \Delta x} \left(x-x_i\right)\hspace{1cm} x_i\leq x\leq x_{i+1}.

Then

\int_{x_i}^{x_{i+1}} f(x) dx =  \int_0^{\Delta x} f(x_i+y) dy \approx f_i \Delta x + (f_{i+1}-f_i){\Delta x\over 2}  = {f_i+f_{i+1}\over 2}\Delta x.

This is known as the trapezoidal rule because it corresponds to the area of the trapezoid formed by connecting the points (x_i,f_i) and (x_{i+1}, f_{i+1}) by a straight line.

When we sum over all intervals, each point gets counted twice except for the left and right boundaries, so

I \approx \Delta x \left[{f_1\over 2} + f_2 + f_3 \dots + f_{N-1} + {f_N\over 2}\right].

Simpson's rule

For the next order, we consider the double interval x_{i-1}<x<x_{i+1}, and write

f(x) \approx f_i + (x-x_i) \left.{df\over dx}\right|_{x_i}+ {(x-x_i)^2\over 2}\left.{d^2f\over dx^2}\right|_{x_i}.

The reason for considering the double interval from x_{i-1} to x_{i+1} is that the linear term is odd in this interval (antisymmetric about x_i) and so does not contribute to the integral, i.e.

\int_{x_{i-1}}^{x_{i+1}} f(x) dx \approx 2 f_i \Delta x + {1\over 2} \left.{d^2f\over dx^2}\right|_{x_i} {2(\Delta x)^3\over 3}.

Taking

\left.{d^2f\over dx^2}\right|_{x_i} \approx {f_{i+1}-2f_i + f_{i-1}\over (\Delta x)^2}

and simplifying then gives

\int_{x_{i-1}}^{x_{i+1}} f(x) dx \approx {\Delta x\over 3} \left[f_{i-1} + 4f_i + f_{i+1} \right].

Adding up the contributions from each double interval, the total integral is

I \approx {\Delta x\over 3}\left[ f_1 + 4f_2 + 2f_3 + 4f_4 + 2f_5 \dots 4f_{N-1} + f_N \right],

where we need the total number of points N to be an odd number (so we can divide the domain into a set of double intervals).

Implement these three methods to calculate the integral 

$$\int_0^{\pi/2} \cos x \, dx = 1.$$

How does the error in each method scale with the number of points $N$? Is it what you expected?

Next, use your code to check the result that the average value of $\sin^2(x)$ is $1/2$, i.e.

$${1\over \pi} \int_{0}^{\pi} \sin^2(x) dx = {1\over 2}.$$

Does the result surprise you? What happened?

What kind of functions are integrated exactly (to machine precision) for each method? (Hint: polynomials of a certain order, which order?)

Gaussian quadrature

In the previous examples, we were given the function f(x) evaluated at a pre-determined set of points \{ x_i\}. But if instead we are able to choose the values of x where we can evaluate the function, we can do better using the method of Gaussian quadratures.

The simplest example is the integral

\int^{1}_{-1} f(x) dx

(which is quite general because for different integration limits you can make a change of variables to bring the limits to -1 and +1). We write the integral as

\int^{1}_{-1} f(x) dx = \sum_{i=1}^N w_i f(x_i).

For a suitable choice of the weights w_i and the locations x_i it can be shown that this expression is exact when f(x) is a polynomial of order 2N-1 or less.

This means that with 2 evaluations of $f(x)$ we can exactly evaluate the integral of a cubic polynomial! How does this compare with Simpson's method?

The proof that this works and the calculation of the values of w_i and x_i is quite involved. [It uses the theory of orthogonal polynomials, so for example for this particular form of the integral, the values of x_i are the roots of the Legendre Polynomial P_N(x).] However, we can look up w_i and x_i using numpy.polynomial.legendre.leggauss

For example,

>>> np.polynomial.legendre.leggauss(2)
(array([-0.57735027,  0.57735027]), array([1., 1.]))

gives the N=2 values. In this case, the weights are both 1, and the x values are \pm 1/\sqrt{3}.

Using [`numpy.polynomial.legendre.leggauss`](https://numpy.org/doc/stable/reference/generated/numpy.polynomial.legendre.leggauss.html)
to get the locations and weights for different choices of $N$, implement Gaussian quadratures. 

Check that the answer is indeed exact for polynomials of degree $2N-1$ or less. What happens for higher order polynomials? 

Try a function that is not a simple polynomial, e.g. $e^{-x^2}$. What is the error in the approximation and how does it scale with $N$?

Hints: 
- to get the correct value of the integral to compare with to see how accurate your approximation is, you could use the general purpose integrator [`scipy.integrate.quad`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html)
- to generate a polynomial with degree `N` and random coefficients between -10 and +10 (for example) you can use
`np.polynomial.Polynomial(np.random.randint(-10,high=10,size=N+1))`

Gaussian quadratures can be applied more generally to integrals of the form

\int W(x) f(x) dx

for some weight function W(x) (in the previous example we had W(x)=1). We write the integral as a sum

\int W(x) f(x) dx= \sum_{i=1}^N w_i f(x_i)

and find the choices of w_i and x_i that make the sum exact for a polynomial of degree 2N-1 or less.

One example is W(x)=e^{-x^2}, ie. integrals

\int_0^\infty e^{-x^2} f(x) dx.

In this case, the weights and locations are given by numpy.polynomial.hermite.hermgauss. [The locations x_i are the roots of the $N$th Hermite polynomial, which are the polynomials that are orthogonal under an inner product defined with weight function W(x).]

Modify your code to use the Gauss-Hermite coefficients and check that you can get an exact answer for the integral of $e^{-x^2}$.

Hint: If you want to use `scipy.integrate.quad` again to get the value of the integral as a comparison, note that you can give it limits of $-\infty$ to $+\infty$ using `-np.inf` and `np.inf`.

Other examples are

• W(x)=e^{-x} with integration limits 0 to \infty. In this case, we need Gauss-Laguerre integration -- see numpy.polynomial.laguerre.laggauss
• W(x)=1/\sqrt{1-x^2} from x=-1 to 1 is Gauss-Chebyshev quadrature -- see numpy.polynomial.chebyshev.chebgauss

Integration challenge

Use Simpson's rule, Gaussian quadrature, and the general purpose integrator [`scipy.integrate.quad`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html) to evaluate the average velocity $\langle\left|v\right|\rangle$ for the 3D [Maxwell-Boltzmann distribution](https://en.wikipedia.org/wiki/Maxwell–Boltzmann_distribution).

For each method, check the numerical error comparing to the analytic result. How many points do you need to get to $0.1$% accuracy?

For Simpson's rule you can use your own implementation from above or you could try [`scipy.integrate.simpson`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.simpson.html).

For Gaussian quadrature, try both Gauss-Hermite and Gauss-Laguerre. Which one is best?

Further reading

• Integration methods are covered in Chapter 7 of Gezerlis.

• Overview of `scipy.integrate'

• QUADPACK is the Fortran 77 library that is used by scipy.integrate.quad under the hood.