Numerical Integration

This notebook follows TAK § 3.3.

Programming concepts introduced:

  • lambda functions.

The goal:

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.

In [1]:
# Standard imports
import numpy as np
import matplotlib.pyplot as plt
import math as m
from mpmath import mp, iv

Lambda functions in Python

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:

In [2]:
lambda x: m.sin(3*x**2)
Out[2]:
<function __main__.<lambda>>

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.

In [3]:
f = lambda x: m.sin(3*x**2)
f(0)
Out[3]:
0.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:

In [4]:
g = lambda x,y,z: (x+y)/z
In [5]:
g(1,2,3)
Out[5]:
1.0

Lambda functions are described in the Python 3 tutorial.

Well-known Quadrature Rules:

The Midpoint Rule:

$$\int_a^b f(x)~dx \approx (b-a)~~f\left(\frac{a+b}{2}\right).$$

The Trapezoid Rule:

$$\int_a^b f(x)~dx \approx \frac{b-a}{2} \big[\,f(a)+f(b)\big].$$

Simpson's Rule:

$$\int_a^b f(x)~dx \approx \frac{b-a}{6} \left[\,f(a)+4 f\left(\frac{a+b}{2}\right) + f(b) \right].$$

Example:

Use Simpson's rule to approximate $$\int_{-\frac{\pi}{2}}^\frac{\pi}{2} \cos x~dx.$$ determine the absolute error in the approximation.

In [6]:
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
Out[6]:
2.0943951023931953
Using calculus we see that $$\int_{-\frac{\pi}{2}}^\frac{\pi}{2} \cos x~dx= [\sin x]_{-\frac{\pi}{2}}^\frac{\pi}{2} = \sin(\frac{\pi}{2})-\sin(\frac{-\pi}{2})=1-(-1)=2.$$ The *absolute error* is the absolute value of the difference between an approximation and the actual value, so:
In [7]:
absolute_error = abs(approx - 2)
absolute_error
Out[7]:
0.09439510239319526

Degree of precision

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}$.

Remark: To check exactness of all polynomials of degree up to $m$, it is enough to check the polynomials $$1, x, x^2, \ldots, x^m.$$ 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:

  • If $f$ is a function and $s \in {\mathbb R}$ is a scalar, then $A(s f)=s A(f)$ where $sf$ denotes the function sending $x$ to $s f(x)$.
  • If $f$ and $g$ are functions, then $A(f+g)=A(f)+A(g)$, where $f+g$ denotes the function sending $x$ to $f(x)+g(x)$. The integral is also linear in the sense above. Using the two rules above, you can show (for example) that if $$A(1) = \int_a^b 1~dx \quad \text{and} \quad A(x) = \int_a^b x~dx,$$ then $A$ is exact for any linear function, i.e., $A(mx+c)=\int_a^b mx+c~dx$ for all $m,c \in {\mathbb R}.$
Result. The Midpoint Rule and Trapezoid Rule have degree of precision $1$.

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.

Result. The Simpson's Rule has degree of precision $3$.

This is explained in TAK on page 54, so we will not discuss it here.

Problem

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:

In [8]:
def A(f):
    return 5/4*f(1) + 15/4*f(3)

Here we expect to get $\int_0^5 1~dx=5$:

In [9]:
A(lambda x: 1)
Out[9]:
5.0

Here we expect to get $\int_0^5 x~dx=12\frac{1}{2}$:

In [10]:
A(lambda x: x)
Out[10]:
12.5

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:

In [11]:
A(lambda x: x**2)
Out[11]:
35.0

Gaussian Quadrature

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).

Changes of coordinates

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:

Change of coordinates formula for quadrature rules. Suppose $A(f)=\sum_{k=0}^N c_k f(x_k)$ is an quadrature rule with degree of precision $m$ over $[-1,1]$. Let $[c,d]$ be another interval and define $g(x)=\frac{d-c}{2} x + \frac{c+d}{2}.$ Let $y_k=g(x_k)$ and $d_k = \frac{d-c}{2} c_k$. Then the quadrature rule $B(f)=\sum_{k=0}^N d_k f(y_k)$ has degree of precision $m$ over $[c,d]$.

We apply this below to the interval from $[0,\pi]$ to estimate an integral.

In [12]:
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$.

In [13]:
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))
For x^0, our approximation is 3.141592653589793 and the exact value is 3.141592653589793.
The error was 0.0
For x^1, our approximation is 4.934802200544679 and the exact value is 4.934802200544679.
The error was 0.0
For x^2, our approximation is 10.335425560099939 and the exact value is 10.335425560099939.
The error was 0.0
For x^3, our approximation is 24.352272758500604 and the exact value is 24.352272758500604.
The error was 0.0
For x^4, our approximation is 61.20393695705627 and the exact value is 61.203936957056285.
The error was -1.4210854715202004e-14
For x^5, our approximation is 160.23153226255067 and the exact value is 160.2315322625507.
The error was -2.842170943040401e-14
For x^6, our approximation is 430.39178495819266 and the exact value is 431.4704611109702.
The error was -1.0786761527775184

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.