Numerical Integration - Trapezium Method

Numerical integration can provide a solution to the problem of integrating functions for which an analytical solution is not available (either because one does not exist or not having enough paper and pencils to derive it).  The trapezium method is the simplest to understand and implement.

Integration is often described as the calculation of the area under the curve.  It the area under the curve is divided into a number of strips, each strip can be approximated to a trapezium and the sum of the area of all the trapeziums will be an approximation of the area under the curve.  As the number of strips increases (and the width of each strip decreases), the approximation will improve.

The algorithm described below, starts with a single strip and then doubles the number of strips until the required level of accuracy is achieved.

Caveat

As with all numerical methods it is important to determine that the algorithm is fit-for-purpose.  In this case the function to be integrated should not be discontinuous or take the value of infinity with the interval over which it is being integrated.  In addition computational efficiency, accuracy and convergence should be considered.

The Area of a Trapezium

The diagram shows a trapezium and the formula for calculating its area.

The formula can be interpreted as the area being the product of the width and the average height.

Algorithm

The algorithm starts with a single strip (Number of strips = 2⁰), then on the next loop the number of strips is increased to 2 (number of strips = 2¹), then four strips (number of strips = 2²) and so on.  On each loop, the approximation to the integral is refined by using the sum of the values of additional ordinates.  After a set number of loops, a check is made on the accuracy.  If the required accuracy is not achieved, the algorithm is deemed to have failed.

Initial Estimate (n=0)

The initial estimate is simply based the start and end points of the interval

The value of the sum of the ordinate heights is expressed in terms of the Integral approximation.

First Loop (n=1)

Next, we have a series of loops, at the start of each loop the number of strips is doubled, on the first loop, this results in a single additional ordinate value being included in the approximation.

Only the value of the new ordinate is calculated, the results of initial estimate are used to provide an updated estimate.

Second Loop (n=2)

The number of strips increases to four

As before, the new ordinate values are used to update the previous estimate

Subsequent loops

The process established in the initial loops, can be generalized to:

Formula

The generalized equation can be developed into a form which can be used as the basis for a computer program

Example code

There are two versions of the code, the first was created in VB.net Express, this is based on a library which is in regular use, the second is in Python which is a translation of the VB.net code.  It consists of two elements, the function trapintegrate which implements the methodology described above and a sample function which is the probability density function of the standard normal distribution.

The table below was derived from Excel.  With the tolerance set to 0.0001 intrap produces similar values to those of Excel when formatted to 4 d.p.

x	Excel
0.0	0.0000
0.2	0.0793
0.4	0.1554
0.6	0.2257
0.8	0.2881
1.0	0.3413
1.2	0.3849
1.4	0.4192
1.6	0.4452
1.8	0.4641
2.0	0.4772
2.2	0.4861
2.4	0.4918
2.6	0.4953
2.8	0.4974
3.0	0.4987
3.2	0.4993
3.4	0.4997
3.6	0.4998
3.8	0.4999
4.0	0.5000
4.2	0.5000
4.4	0.5000
4.6	0.5000
4.8	0.5000

Reducing the value of tolerance, imporves the accuracy at the expense of increasing the run time.

VB.Net

Python

Python code blocks are defined by indententation, some adjustment might be needed after a copy and paste operation.