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

VB.Net code

The sample code consists of three elements, a call to the integration function which provides the integration interval, a sample function which is the probability density function of the standard normal distribution and the integration function itself which is called TrapIntegrate

Call to TrapIntegrate

'This sample integrates the probability distribution function

'of the standard normal distribution over the range 0 to 1,

'the value of which to 4dp is 0.3413.

Dim strText As String

Dim dblA As Double = 0

Dim dblB As Double = 1

Dim dblI As Double

If TrapIntegrate(dblA, dblB, dblI) Then

strText = String.Format("{0:0.0000}", dblI)

Else

strText = "TrapIntegrate returned false"

End If

Sample Function (called within TrapIntegrate)

Private Function f(ByVal dblX As Double) As Double

'Test function (probability density of standard normal distribution)

Return Math.Exp(-dblX * dblX / 2) / Math.Sqrt(2 * Math.PI)

End Function

Integration Function

Private Function TrapIntegrate(ByVal dblA As Double, _
                                            ByVal dblB As Double, _
                                            ByRef dblI As Double) As Boolean

'TrapIntegrate uses the integrates a function over the interval

'A to B using the trapezoidal method. Returns true if the

'required level of accuracy is achieved, otherwise false. The

'variable I contains the value of the integral

Dim intN As Integer = 1 '2^N Number of strips

Dim intJ As Integer 'Loop counter

Dim intK As Integer 'Convenience variable - 2^(N-1)

Dim intNMin As Integer = 5 'Min number of strips (2^Nmin)

Dim intNMax As Integer = 10 'Max number of strips (2^Nmax)

Dim dblTol As Double = 0.0001 'Fractional error

Dim dblFSum As Double 'Running total of function values

Dim dblIPrev As Double 'Integral from previous loop

'Initial estimate of integral

dblIPrev = (dblB - dblA) * (f(dblA) + f(dblB)) / 2

'Continuously double the number of strips until required

'accuracy is achieved or number of strips becomes

'stupidly large

Do While intN <= intNMax

'Reset Running total of function values

dblFSum = 0

'Convenience variables

intK = Math.Pow(2, intN - 1)

'sum of additional function values

For intJ = 1 To intK

dblFSum = dblFSum + f(dblA + (dblB - dblA) * (2 * intJ - 1) / 2 / intK)

Next

'Revise value of Integral

dblI = 0.5 * (dblIPrev + (dblB - dblA) * dblFSum / intK)

'Check accuracy

If intN >= intNMin Then

If Math.Abs(dblI - dblIPrev) < dblTol * Math.Abs(dblIPrev) Or _
  (dblI = 0 And dblIPrev = 0) Then

Return True

End If

End If

'Preserve value of integral for testing on next loop

dblIPrev = dblI

'Increment N (Double the number of strips)

intN = intN + 1

Loop

'No convergence

Return False

End Function

Page updated: 26-Jun-08

For more information: info@brighton-webs.co.uk