This notebook is a modified and annotated version of the example shown in class November 26, 2015.

It demonstrates how to implement the trapezoidal rule for computing an approximate value for a definite integral.

In [1]:
%pylab

Using matplotlib backend: TkAgg
Populating the interactive namespace from numpy and matplotlib

We set up a simple example: Find the integral of the function f(x)=x^2 on the interval [0,1].

To make the prodecure more general (cf. Homework Set 8), we take integration nodes which are not equidistantly spaced. The following cell just sets up random nodes in the interval [0,1]. The details of this do not matter for the purpose of understanding the trapezoidal rule and shall not be explained further.

In [2]:
x=cumsum(rand(10))
x=r_[0,x]/x[-1]
x
Out[2]:
array([ 0.        ,  0.10197122,  0.24503111,  0.30551511,  0.33595835,
        0.37408344,  0.4849629 ,  0.5901724 ,  0.69152883,  0.86217171,  1.        ])

Now let's compute an array containing the corresponding values of the function f(x)=x^2:

In [3]:
f=x**2
f
Out[3]:
array([ 0.        ,  0.01039813,  0.06004024,  0.09333948,  0.11286801,
        0.13993842,  0.23518902,  0.34830347,  0.47821213,  0.74334006,  1.        ])

For the trapezoidal rule, we need to compute the width of each integration bin. An array containing the upper ends of the integration bins is provided by

In [4]:
x[1:]
Out[4]:
array([ 0.10197122,  0.24503111,  0.30551511,  0.33595835,  0.37408344,
        0.4849629 ,  0.5901724 ,  0.69152883,  0.86217171,  1.        ])

An array containing the corresponding lower ends of the integration bins is provided by

In [5]:
x[:-1]
Out[5]:
array([ 0.        ,  0.10197122,  0.24503111,  0.30551511,  0.33595835,
        0.37408344,  0.4849629 ,  0.5901724 ,  0.69152883,  0.86217171])

Therefore, an array with the bin widths is computed as

In [6]:
dx=x[1:]-x[:-1]

In exactly the same way, an array with the average height of the function in each bin is given by

In [7]:
h=(f[1:]+f[:-1])/2

The array containing the area of each bin is therefore obtained by

In [8]:
dx*h
Out[8]:
array([ 0.00053015,  0.00503845,  0.00463851,  0.00313881,  0.00481913,
        0.02079696,  0.03069448,  0.04188634,  0.10422459,  0.12014079])

The total area, approximating the integral, is therefore given by

In [9]:
sum(dx*h)
Out[9]:
0.33590821838797308

The Fundamental Theorem of Calculus tells us that the exact area under the graph of f on [0,1] is

In [10]:
1/3.
Out[10]:
0.3333333333333333

We see that the computed value is rather close, but not exact. The approximation will get better by taking more bins, i.e, replacing the size parameter 10 in Cell 2 by a larger value.