Trapezoidal Rule

The Trapezoidal rule is a mathematical method for doing numerical integration ("calculating the area under a curve"). This is often useful when an exact integral does not exist, can not easily be obtained, or is mathematically too time consuming for repetitious automated calculations.

One approach to obtain a numerical solution of an integral is to approximate the function with an n^th order polynomial, since these are relatively simple to integrate. The choice of the order of the polynomial depends on the required accuracy (higher order generally results in a higher precision), and the number data points over the selected interval (higher order requires more data points).

Consider the function f(x). We want to calculate the area under the curve over the interval a ≤ x ≤ b:

         b
     I = ∫ y(x)dx
         a

The simplest polynomial approximation, a 0^th order polynomial would draw a straight horizontal line, halfway between y(a) and y(b), marking a rectangle. The area of this rectangle is simply the length times the height. But that is not a very accurate method; the approximation of the area yields a relatively large error.

A better approximation is to approximate the curve with a 1^st order polynomial; a line that goes trough the points a and b. This can be seen in the following figure:

   y(x)

y(b)_|                  ..
     |               ./|C
     |              ./ |
     |             ./  |
     |             /   |           .....  y(x)
     |            /.   |
     |           / .   |
     |          / .    |
     |       . /.      |
y(a)_|     . D|        |
     |    .   |        |
     | ..     |        |
     |.       |        |
     |        |        |
     |_______A|________|B___ x
              a        b

Over the range a ≤ x ≤ b, the function y(x) is now approximated with a trapezoid marked by ABCD. The area of the trapezoid is:

     ABCD = 1/2 AB(AD + BC) = 1/2 (b-a)(y(a) + y(b))

We can gain a higher precision if we divide the interval a ≤ x ≤ b up into several equal parts, a = x₀, x₁, x₂, ..., x_n-1, x_n = b. We find the corresponding y values, y₀ = y(x₀), y₁ = y(x₁), y₂ = y(x₂), ..., y_n = y(x_n). Each part has a width equal to h = (b - a)/n

The Trapezoidal Rule for approximating the integral now becomes:

     b
     ∫ y(x)dx = h/2(y0 + 2y1 + 2y2 + ... + 2yn-1 + yn)
     a

Example: Let's approximate the integral of the function:

     y = x2 + 2x -3

over the range 3 ≤ x ≤ 4. The exact solution is straightforward:

     4                            4
     ∫ y(x)dx = [ 1/3 x3 + x2 - 3x] = 76/3 - 9 = 16.3333...
     3                            3

Now let's try the numerical solution. First generate a table with x-, and y-values for the selected interval:

 n       xn      y
 0      3.0    12.00   
 1      3.2    13.64
 2      3.4    15.36
 3      3.6    17.16
 4      3.8    19.04
 5      4.0    21.00

The numerical solution is given by (h=1/5):

     4                
     ∫ y(x)dx = (1/10)(12.00 + 2×13.64 + 2×15.36 + 2×17.16 + 2×19.04 + 21.00) = 16.34
     3                

The example shows that the trapezoidal rule can obtain relatively accurate results (in this case an error of approximately 0.04%) with only a small number of simple calculations. The error between the exact integral (I) and the numerical solution I_n is given by:

     I - In = -1/12 h3y"(ε)

with a ≤ ε ≤ b

Of course, the integral can also be calculated with higher order polynomial approximations. The 2^nd order polynomial approximation is called Simpson's Rule.