This notebook follows TAK § 3.3.
Programming concepts introduced:
Approximate an integral $$I = \int_a^b f(x)~dx$$ by a weighted sum of function values. That is, we choose points $x_0, \ldots, x_N$ in the interval $[a,b]$ and weights $c_0, \ldots, c_N \in {\mathbb R}$ and hope that $$I \approx \sum_{k=0}^N c_k f(x_k).$$ (Here $\approx$ denotes approximately.) The sum on the right is called a quadrature rule.
# Standard imports
import numpy as np
import matplotlib.pyplot as plt
import math as m
from mpmath import mp, iv
In Python we can create functions in two ways. We have covered the first way. The second ways makes it easy to create functions in one line. It also allows you to define a function without giving it a name.
We have this notion in mathematics two. For example, $$x \mapsto \sin(3 x^2)$$ defines a function without giving it a name. The domain is implicit. This function can be defined in python using the notation:
lambda x: m.sin(3*x**2)
You can see that the object returned by the statement above is a function object. In order to use it in the usual way, we should give it a name.
f = lambda x: m.sin(3*x**2)
f(0)
The line f = lambda x: m.sin(3*x**2) is equivalent to the code:
def f(x):
return m.sin(3*x**2)We can also define functions which take more than one variable as input:
g = lambda x,y,z: (x+y)/z
g(1,2,3)
Lambda functions are described in the Python 3 tutorial.
Use Simpson's rule to approximate $$\int_{-\frac{\pi}{2}}^\frac{\pi}{2} \cos x~dx.$$ determine the absolute error in the approximation.
def f(x):
return m.cos(x)
a = -m.pi/2
b = m.pi/2
approx = (b-a)/6 * (f(a) + 4*f((a+b)/2) + f(b))
approx
absolute_error = abs(approx - 2)
absolute_error
An integration rule over an interval $[a,b]$ has degree of precision $m$ if it is exact (zero error) for all polynomials of degree at most $m$, but is not exact for $x^{m+1}$.
This is because of our choice of approximation as a weighted sum, $$A(f) = \sum_{k=0}^N c_k f(x_k).$$ This expression is linear in $f$, meaning that:
Proof for the midpoint rule: Let $A$ denote the midpoint rule applied to the interval $[a,b]$. We need to show that $A$ is exact for the function $1$, and the function $x$, but not for $x^2$.
First we check it for the constant function $1$. Observe $\int_a^b 1~dx=b-a$ and $$A(1)=(b-a)\cdot 1=b-a.$$ So $A$ is exact for $1$.
Now we check it for the function $x$. Observe that $$\int_a^b x~dx=\left[\frac{1}{2}x^2\right]_a^b = \frac{1}{2}(b^2-a^2).$$ We also have $$A(x) = (b-a) \cdot \frac{b-a}{2}=\frac{1}{2}(b^2-a^2).$$ Thus, $A$ is exact for $x$.
To see that $A$ is not exact for $x^2$, observe that $$\int_a^b x^2~dx=\left[\frac{1}{3}x^3\right]_a^b = \frac{1}{3}(b^3-a^3)$$ and $$A(x^2) = (b-a) \cdot \left(\frac{b-a}{2}\right)^2=\frac{1}{4}(b-a)^3.$$ QED
We will not prove the trapezoid case. It should be clear geometrically.
This is explained in TAK on page 54, so we will not discuss it here.
Derive a rule which has degree of precision $1$ using the points $1$ and $3$ in the interval $[0,5]$.
Solution.
A rule using points $1$ and $3$ must have the form $$A(f) = c_1 f(1) + c_3 f(3).$$ Since we want our rule to have degree $1$, it must satisfy $$A(1)=\int_0^5 1~dx =5 \quad \text{and} \quad A(x)=\int_0^5 x~dx =\frac{25}{2}.$$ Observe that by definition, we have $$A(1)=c_1 + c_3 \quad \text{and} \quad A(x) = c_1 + 3 c_3.$$ So, we must solve the equations $$c_1+c_3 = 5 \quad \text{and} \quad c_1 + 3 c_3=\frac{25}{2}.$$ Subtracting the first equation from the second yields $$2 c_3 = \frac{15}{2} \quad \text{so} \quad c_3=\frac{15}{4}.$$ Then plugging this value of $c_3$ into one of the equations and solving yields $$c_1=\frac{5}{4}.$$ Thus, our approximation rule has the form $$A(f) = \frac{5}{4} f(1) + \frac{15}{4} f(3).$$
Below we check that our formula works as expected:
def A(f):
return 5/4*f(1) + 15/4*f(3)
Here we expect to get $\int_0^5 1~dx=5$:
A(lambda x: 1)
Here we expect to get $\int_0^5 x~dx=12\frac{1}{2}$:
A(lambda x: x)
We can see that our formula is not exact in degree $2$ since $\int_0^5 x^2~dx=\frac{125}{3}=41\frac{2}{3}$, but:
A(lambda x: x**2)
The idea of Gaussian Quadrature is to find the approximation with the highest degree of precision among all approximations using $N$ points. It is know that you can always find an approximation with degree $2N-1$.
The book works out the formula for the Gaussian Quadrature using $2$ and $3$ points in the interval $[-1,1]$. The two point formula is $$G_2(f) = f\left(\frac{1}{\sqrt{3}}\right) + f\left(\frac{-1}{\sqrt{3}}\right),$$ which has degree of precision $3$. The three-point formula is $$G_3(f) = \frac{5}{9}\,f\left(\sqrt{\frac{3}{5}}\right) + \frac{8}{9}\,f(0) + \frac{5}{9}\,f\left(-\sqrt{\frac{3}{5}}\right),$$ which has degree of precision $5$. These are derived in example 9 (starting on page 55 of TAK).
Recall the change of coordinates formula for integration. If $f$ is integrable and $g$ is differentiable, then $$\int_a^b f\big(g(x)\big)g'(x)~dx = \int_{g(a)}^{g(b)} f(x)~dx.$$ We will apply this to the $3$-point Gaussian Quadrature to derive a formula which holds over any interval.
We can apply this to move an integral defined over say $[-1,1]$ to any other interval. If we want to consider an integral defined over $[c,d]$ we can take $$g(x)=\frac{d-c}{2} x + \frac{c+d}{2}.$$ This is a linear function, thus if $f(x)$ is a polynomial $f\big(g(x)\big)g'(x)$ is a polynomial of the same degree. Applying our change of coordinates formula, and pulling out the constant $g'(x)=\frac{d-c}{2}$ yields $$\frac{d-c}{2} \int_{-1}^1 f\big(g(x)\big)~dx = \int_{c}^{d} f(x)~dx.$$ The integral on the left is over the interval $[-1,1]$ so our $G_3$ will be exact provided $f\big(g(x)\big)$ is a polynomial of degree 5 or less. Thus for polynomials $f$ of degree $5$ or less, we have $$\frac{d-c}{2} \left[\frac{5}{9}\,f \circ g\left(\sqrt{\frac{3}{5}}\right) + \frac{8}{9}\,f \circ g(0) + \frac{5}{9}\,f\circ g\left(-\sqrt{\frac{3}{5}}\right)\right] = \int_{c}^{d} f(x)~dx.$$ We find it easiest to leave our expression like this, but we could expand $g(\star)$ to get points, so this is still in the form of an integration rule. Also it still has degree of precision $5$. This is a special case of a more general rule:
We apply this below to the interval from $[0,\pi]$ to estimate an integral.
c = 0
d = m.pi
# Linear change of coordinates [-1,1] -> [c,d]
g = lambda x: (d-c)/2*x + (c+d)/2
# Points used in Gaussian Quadrature over [-1,1]
x1 = -m.sqrt(3)/m.sqrt(5)
x2 = 0
x3 = m.sqrt(3)/m.sqrt(5)
# Weights used in Gaussian Quadrature over [-1,1]
c1 = 5/9
c2 = 8/9
c3 = 5/9
# Points used in Gaussian Quadrature over [c,d]
y1 = g(x1)
y2 = g(x2)
y3 = g(x3)
# Weights used in Gaussian Quadrature over [c,d]
d1 = (d-c)/2 * c1
d2 = (d-c)/2 * c2
d3 = (d-c)/2 * c3
# Gaussian Quadrature formula over [c,d]
A = lambda f: d1*f(y1) + d2*f(y2) + d3*f(y3)
Our $A$ should be have degree of precision $5$. We will check that.
Below we compare $A(x^k)$ and $\int_0^\pi x^k~dx=\frac{1}{k+1}\pi^{k+1}$ for values of $k \leq 6$.
for k in range(7):
approx = A(lambda x: x**k)
exact = m.pi**(k+1) / (k+1)
print("For x^{}, our approximation is {} and the exact value is {}.".format(k, approx, exact))
print("The error was {}".format(approx-exact))
It seems clear that it was exact up to degree $5$ (with non-zero values of the order $10^{-14}$ probably coming from errors in floating point arithmetic.) The degree $6$ error of about $1$ indicates that our formula is not exact for polynomials of degree $6$, which we expect.