Precept 11: Numerical integration

We experiment with numerical quadrature to estimate:

$$ \int_{-1}^1 f(x) dx $$

Implement the Trapezoidal rule:

$$ T(h) = - \frac{h}{2} ( f(x_0) + f(x_n)) + h\sum_{i=0}^{n} f(x_i) $$

where $x_i = -1 + \frac{2i}{n}$ , $h = \frac{2}{n}$, $i = 0, \dots n$.

import numpy as np
def trapezoidal_rule(f,n):
    return 0

1. Test your rule on a smooth function:

$$g(x) = 1/ (1 + 25 x^2) $$

Make a log-log plot displaying the $O(h^2)$ convergence

import matplotlib.pyplot as plt
g = lambda x : 1 / ( 1 + 25 * x**2)

analytical_value = 2/5 *  np.arctan(5)
nns = np.logspace(1,4, 20).astype(int)
errors = np.zeros(nns.size)
for n_idx, n in enumerate(nns):
    a = 1
    # Add your code here
plt.loglog( )

[]

1. Test it on a smooth and periodic function:

$$ g(x) = \exp( \sin(2 \pi x )) $$

Make a semilog plot illustrating the spectral accuracy.

import matplotlib.pyplot as plt
g2 = lambda x : np.exp(np.cos(2*np.pi*x))
analytical_value = 2.5321317555040166711964892504294350752153406227099244136162706624
nns = np.linspace(1,40).astype(int)
errors = np.zeros(nns.size)

## Add your code here

plt.semilogy()

[]

1. Use the Richardson extrapolation to get a $O(h^4)$ quadrature scheme based on the Trapezoidal rule:

$$T_2(h) = \frac{4 T(h/2) - T(h)}{3}$$

Illustrate the convergence rate on the function

$$g(x) = 1/ (1 + 25 x^2) $$

a log-log plot.

import matplotlib.pyplot as plt
g = lambda x : 1 / ( 1 + 25 * x**2)
analytical_value = 2/5 *  np.arctan(5)
nns = np.logspace(1,4, 20).astype(int)
errors = np.zeros(nns.size)
## add your code here


plt.loglog()

[]

1. Legendre-Gauss quadrature.

Let $x_1 \le \cdots \le x_n$ be the roots of the $n$-th degree Legendre Polynomial $P_n(x)$, and let $w_1,\ldots,w_n$ be weights such that $$ \sum_{k=1}^n P_j(x_k) w_k = \int_{-1}^1 P_j(x) dx $$ for $j = 0,\ldots,n-1$. Together, the points $x_1,\ldots,x_n$ and weights $w_1,\ldots,w_n$ form a quadrature rule, which is a numerical rule for integrating a function. Given a function $f : [-1,1] \rightarrow \mathbb{R}$ we hope to approximate the integral of $f$ over $[-1,1]$ by $$ \int_{-1}^1 f(x) dx \approx \sum_{k=1}^n f(x_k) w_k $$

The weights and nodes can be computed by: [x,w] = np.polynomial.legendre.leggauss(n)

This quadrature is spectrally accurate for smooth function. Make a semilog plot showing this behavior.

import matplotlib.pyplot as plt
g = lambda x : 1 / ( 1 + 25 * x**2)
analytical_value = 2/5 *  np.arctan(5)
nns = np.linspace(1, 100).astype(int)
errors = np.zeros(nns.size)

## add your code here

plt.semilogy()

[]

1. The degree-n Gauss-Legendre quadrature is exact for polynomials of degree $2n-1$.

Check that your function exactly integrates $ f(x) = (x^2+1)(x-3)(x-1)^2$ on $[-1,1]$ when $n=3$.

import matplotlib.pyplot as plt
g = lambda x : (x**2 +1)*(x-3)*((x-1)**2)

analytical_value = - 40/3
## add your code ehre
gl_quad = 0
error = np.abs(analytical_value - gl_quad)
print(error)

13.333333333333334