Numerical Integration

We are going to look at three different strategies for approximating a definite integral:

  • The Rectangle Method

  • The Trapezium Method

  • Simpson's Method

All these methods are concerned with approximating the integral:

$$ \int_a^b f(x) dx\,, $$
where $f(x)$ is a given function and $a$ and $b$ are given numbers.

Recall that this integral represents a certain area, as seen in the sketch included below





The definite integral above represents the area between the graph of the function and the x-axis, which is bounded by the vertical lines x=a and x=b

The rectangle Method

This is based on approximating the aforementioned area by using rectangles, as seen below.

We are in control regarding how many rectangles we want to use. It should be intuitively clear that our area is going to be well approximated by having a very large number of such rectangles.

For simplicity, we shall consider a set of equidistant points between $x=a$ and $x=b$:

$$ a = x_0 < x_1 < x_2 <\dots < x_{n-1} < x_n = b\,. $$
'Equidistant' means that $x_{i+1}-x_i = h$ (constant) for all $i=0,1,\dots, n-1$.





The Rectangle Method

Let's call the various rectangles: $R_0$, $R_1,\,\dots,,R_{n-1}$, and their corresponding areas $I_0$, $I_1,\,\dots,,I_{n-1}$. Then

\begin{alignat*}{1} &I_ 0 = (x_1-x_0)f(x_0) = hf(x_0)\,,\\[0.2cm] &I_1 = (x_2-x_1)f(x_1) = hf(x_1)\,,\\[0.2cm] &\dots\,,\\[0.2cm] &I_k = (x_{k+1} - x_k)f(x_k) = hf(x_k)\,,\\[0.2cm] &\dots\,,\\[0.2cm] &I_{n-1} = (x_n-x_{n-1})f(x_{n-1}) = hf(x_{n-1}) \end{alignat*}
and hence

\begin{equation*} \large{\int_a^bf(x)\, dx \simeq I_0 + I_1+ I_2+\dots + I_{n-1} =h\sum_{k=0}^{n-1} f(x_k)} \end{equation*}
The absolute error for this approximation is

$$ \large A.E. = \mathcal{O}\left(\frac{(b-a)^2}{2n}\max_{a\leq x\leq b} |f'(x)|\right) $$

In [2]:
def integrateREC(fname, xmin, xmax, intervals):
    """ Integrate by using the rectangle rule."""
    h = (xmax - xmin) / float(intervals)
    total = 0.0

    # Perform the integration:
    x = xmin
    for interval in range(intervals):
        # Add the area in the trapeziod for this slice:
        total += h * fname(x)

        # Move to the next slice:
        x += h

    return total  # this returns the approx. of the integral

EXAMPLE:

a) $$ \int_0^{\pi/2} \sin (401 x)\,dx = \frac{1}{401} \simeq 0.00249376 $$

In [19]:
import numpy as np

def f(x): 
    return np.sin(401.0*x)

intervals = [1000, 5000, 10000, 100000]
print("The integral is approximately equal to:")
print(" ")
for interval in intervals:
    print(integrateREC(f, 0, 0.5*np.pi, interval), "\t (n={})".format(interval))
    

The integral is approximately equal to:
 
0.0016253646698 	 (n=1000)
0.00233338698793 	 (n=5000)
0.00241440119218 	 (n=10000)
0.00248590335823 	 (n=100000)

It should be clear from the little experiment above that the rectangle method is terrible....

The Trapezium Method

In this case we divide the interval $[a, b]$ into a set of $n$ intervals, each of length $h$. The function whose definite integral we aim to find is approximated by a straight line in each sub-interval, as shown in the sketch below. Hence the area that represents our integral will be approximately equal to the sum of the areas of the trapeziums seen in that sketch.





The Trapezium Method

Let's call the various trapeziums: $T_0$, $T_1,\,\dots,,T_{n-1}$, and their corresponding areas $I_0$, $I_1,\,\dots,,I_{n-1}$. Then

\begin{alignat*}{1} &I_ 0 = \frac{1}{2}h[f(x_1)+f(x_0)] \,,\\[0.2cm] &I_1 = \frac{1}{2}h[f(x_1)+f(x_2)]\,,\\[0.2cm] &\dots\,,\\[0.2cm] &I_k = \frac{1}{2}h[f(x_k)+f(x_{k+1})]\,,\\[0.2cm] &\dots\,,\\[0.2cm] &I_{n-1} = \frac{1}{2}h[f(x_{n-1})+f(x_n)] \end{alignat*}
and hence

\begin{alignat*}{1} \large{\int_a^bf(x)\, dx }&\large{\simeq I_0 + I_1+ I_2+\dots + I_{n-1}}\\[0.1cm] &\large =\frac{1}{2}h[f(x_0) + 2f(x_1) + 2f(x_2) +\dots + 2f(x_{n-1}) + f(x_n)]\\[0.1cm] &\large=\frac{h}{2}[f(x_0) + f(x_n)] + h\sum_{k=1}^{n-1}f(x_k) \end{alignat*} The absolute error for this approximation is

$$ \large A.E. = \mathcal{O}\left(\frac{(b-a)^3}{12n^2}\max_{a\leq x\leq b} |f''(x)|\right) $$

In [6]:
def integrateTRAPZM(fname, xmin, xmax, intervals):
    """ Integrate by using the trapezoid rule."""
    h = (xmax - xmin) / float(intervals)
    total = 0.0

    # Perform the integration:
    x = xmin
    for interval in range(intervals):
        # Add the area in the trapeziod for this slice:
        total += h * (fname(x) + fname(x + h))/float(2.0)

        # Move to the next slice:
        x += h

    return total  # this returns the approx. of the integral

EXAMPLE:

a) $$ \int_0^\pi e^x\cos x\,dx = -\frac{1}{2}(e^\pi+1)\simeq = -12.070346 $$

In [20]:
import numpy as np

def f(x):
    return np.exp(x)*np.cos(x)

intervals = [10, 25, 50, 200, 400, 600]
print("The integral is approximately equal to:")
print(" ")
for interval in intervals:
    print(integrateTRAPZM(f, 0, np.pi, interval), "\t (n={})".format(interval))
    

The integral is approximately equal to:
 
-12.2695456708 	 (n=10)
-12.1021309039 	 (n=25)
-12.0782893309 	 (n=50)
-12.0708426936 	 (n=200)
-12.0704704099 	 (n=400)
-12.070401469 	 (n=600)

Much better than the rectangle method....

Simpson's Rule

Again, we start with a set of equidistant points in the interval $[a, b]$.

The idea is to approximate the function $f(x)$ by a parabola on each sub-interval $[x_0, x_2]$, $[x_2, x_4]$, etc. For this to be possible $n$ must be an EVEN number (any number that is divisible by 2 is an even number).





Simpson's Method

You have the details of the derivation on the slides posted on Brightspace (lots of details can be found on the internet including YouTube and Khan Academy). Here I'll mention only the final result that is needed for coding this approximation:

$$ \large{{\int_{{a}}^{{b}}}{f{{\left({x}\right)}}}{d}{x}\simeq\frac{{h}}{{3}}{\left({y}_{{0}}+{4}{y}_{{1}}+{2}{y}_{{2}}+{4}{y}_{{3}}+{2}{y}_{{4}}+\ldots+{2}{y}_{{{n}-{2}}}+{4}{y}_{{{n}-{1}}}+{y}_{{n}}\right)}} $$
where you may recall that

$$ \large{h = (b-a)/n\,.} $$
Note the pattern of the coefficients: $$ 1, 4, 2, 4, 2, 4, 2, \dots 4, 2, 4, 2, 4, 2, 4, 1 $$
The absolute error for this approximation is

$$ \large A.E. = \mathcal{O}\left(\frac{(b-a)^5}{180n^4}\max_{a\leq x\leq b} |f''''(x)|\right) $$

In [8]:
def integrateSIMP(f, a, b, N=50):
    """
    Approximate the integral of f(x) from a to b by Simpson's rule.

    INPUTS:
    ----------
    f : function
             of a single variable
    a , b : numbers
        Interval of integration: [a,b]
    N : (even) integer
        Number of subintervals of [a,b]

    OUTPUTS:
    -------
    float
        Approximation of the integral of f(x) over [a, b] 
        using Simpson's rule with N subintervals of equal length.
    """
    if N % 2 == 1:
        raise ValueError("N must be an even integer.")
    h = (b-a)/float(N)          # the stepsize 
    x = np.linspace(a, b, N+1)   # define the mesh 
    y = f(x)
    var = (h/3.0) * np.sum(y[0:-1:2] + 4.0*y[1::2] + y[2::2])  # numerical approx. for integral 
    return var

EXAMPLES:

a).
$$ \int_0^\pi\sin(x)\,dx = 2 $$

In [13]:
integrateSIMP(lambda x: np.sin(x), 0.0, np.pi, 30)
Out[13]:
2.0000013379479493



b).
$$ \int_0^{\sqrt{\pi}} 2x^2\cos(x^2)\,dx $$

In [14]:
def f(x): 
    return 2.0*(x**2)*np.cos(x**2)

integrateSIMP(f, 0, np.sqrt(np.pi), 60)
Out[14]:
-0.89482976891285448



c).
$$ \int_0^{\pi/4}\ln(1+\tan x)\,dx $$

In [15]:
def f(x):
    return np.log(1.0+np.tan(x))

integrateSIMP(f, 0, 0.25*np.pi, 80)
Out[15]:
0.27219826128795022

Trapezium Method with control of accuracy

If we want to approximate the value of an integral to within a specified number of significant digits we can keep repeting the process with $h\to h/2$ (which is of course equivalent to doubling the number of 'intervals'). A possible implementation of such a code is included below.

As a typical example, we shall consider the following integral

$$ \int_{\pi/6}^{\pi/2}\sin^2(x)\cos(x)\,dx = \frac{7}{24}\simeq 0.29166666666666\dots\,. $$
(This integral can be calculated very easily by a simple change of variable.)

In [12]:
import numpy as np

# define the integrand; you can also use a lambda statement:
def f(x):
    return (np.sin(x)**2)*np.cos(x)
#
#
# specify the integration range as well as the tolerance
# if you want 3 significant digits, TOL = 10**(-3)
# we want 6 significant digits, so TOL = 10**(-6)
TOL = 1.0e-6
a = np.pi/6.0
b = np.pi/2.0
#
# run the Trapezium integrator for intervals = 2
intervals = 2
I_old = 0.0
I_new = integrateTRAPZM(f, a, b, intervals)
#
# keep comparing the difference between the two approximations:
while abs(I_new - I_old) >= TOL:
    I_old = I_new # save the I_new calculated at the previous step        
    intervals = 2*intervals  # double the number of sub-intervals
    I_new = integrateTRAPZM(f, a, b, intervals) # find the new approximation
#    
# when we break out of the loop I_new will contain 
# the desired result:    
print("Integral =", I_new)
print("Number of intervals =", intervals)

Integral = 0.291666525045
Number of intervals = 1024

The Midpoint Rule

This is a version of the Rectangle Method (more or less). It is based on approximating the area mentioned above by using a special choice of rectangles. In particular, the height of these rectangles corresponds to

$$ f(x_{i}^*)\,, $$
where

$$ x_i^* = \frac{1}{2}(x_i+x_{i+1})\quad\mbox{is the midpoint of the interval}\qquad (x_i, x_{i+1}) $$
for $i=0,1,2\,\dots, (n-1)$.





The Rectangle Method vs. the Midpoint Method

Let's call the various rectangles: $R_0$, $R_1,\,\dots,,R_{n-1}$, and their corresponding areas $I_0$, $I_1,\,\dots,,I_{n-1}$. Then

\begin{alignat*}{1} &I_ 0 = (x_1-x_0)f(x_0^*) = hf\left(\frac{x_0+x_1}{2}\right)\,,\\[0.2cm] &I_1 = (x_2-x_1)f(x_1^*) = hf\left(\frac{x_1+x_2}{2}\right)\,,\\[0.2cm] &\dots\,,\\[0.2cm] &I_k = (x_{k+1} - x_k)f(x_k^*) = hf\left(\frac{x_k+x_{k+1}}{2}\right)\,,\\[0.2cm] &\dots\,,\\[0.2cm] &I_{n-1} = (x_n-x_{n-1})f(x_{n-1}^*) = hf\left(\frac{x_{n-1}+x_n}{2}\right) \end{alignat*}
and hence

\begin{equation*} \large{\int_a^bf(x)\, dx \simeq I_0 + I_1+ I_2+\dots + I_{n-1} =h\sum_{k=0}^{n-1} f(x_k^*)} \end{equation*}
The absolute error for this approximation is

$$ \large A.E. = \mathcal{O}\left(\frac{(b-a)^3}{24n^2}\max_{a\leq x\leq b} |f''(x)|\right) $$



HOMEWORK:

Write a Python code that implements this method. Test your code on the examples of the various integrals that apear above.