Friday, 6 May 2011

Tutorial: Trapezium Method

A subject I am particularly interested in is everything to do with numerical and mathematical algorithms. For this reason I decided that I am going to start a series on numerics. A second series will be on algorithms and discrete maths. The topics I am going to cover will not necessarily be published in any logical order as I will write about stuff that comes to my mind. But I will provide a page here which will act as a kind of index. In this way I hope that, after some time, I will be able to provide a tutorial on the subjects together with the odd little gem strewn in between.

Trapezium Method of Numerical Integration

When evaluating an integral using the trapezium rule the function f to be integrated is evaluated at discrete points xi between the integration boundaries a and b

xi = a + hi

where the interval size h is given by

h = \frac{b-a}{N},

N is the number of intervals and i=0\ldots N. I will also refer to N as the resolution. In practice N is often chosen to be a power of 2. There is no particular reason for this other than the fact that one commonly increases the resolution by a factor of 2 when improving the accuracy.

Let's abbreviate the notation by defining the discrete function values

f^N_i = f(x_i)

Here I used the upper index N to clarify that the discrete function values depend on the resolution. With this definition we can now write the trapezium approximation TN of resolution N to the integral as

\int_a^b f(x) dx \approx \frac{1}{2h}\left( f^N_0 + 2 \sum_{i=1}^{N-1} f^N_i + f^N_N \right) =: T_N

What this means geometrically is explained in the following graph. The integral of a function is the area below the graph of that function. The trapezium rule approximates this area by adding the areas of the trapezoids created by the x-axis and the function values at the discrete points. Each trapezoid has the area

\frac{1}{2h}\left( f_i + f_{i+1} \right)

Summing these areas up, one obtains the above formula. Let's do an example. We want to integrate the function sin(πx) between the limits 0 and 1. Using N=4 we get

T_4 = \frac{1}{8}\left( 0 +  1.41421 + 2.00000 + 1.41421 + 0 \right) =  0.60355

Comparing this with the exact solution

\int_0^1\sin(\pi x) dx = \frac{2}{\pi}

We see that, with this crude approximation, the error is about \varepsilon=0.03307

Of course one expects this approximation to get better as N increases and thus h decreases. In fact, one can show that, for sufficiently smooth functions

\int_a^b f(x) dx = T_N + O(h^2).

The above equation means that the error in the trapezium approximation decreases quadratically with the number of steps.

Let's see this on our example. In the diagram I plotted the error against the resolution on a double logarithmic plot. Up to N = 1e6 the curve is almost a straight line with a slope of -2. To show this the green line is the function \varepsilon \propto N^{-2}. This means that, in order reduce the error by a factor of four, one has to double the resolution.

One interesting feature of the trapezium rule is the fact that it makes the doubling of the resolution easy for us. When going from N to 2N, one does not have to evaluate the function at all the 2N points. Instead one can use the result TN of the previous approximation. The rule is

T_{2N} = \frac{1}{2} T_N + \frac{1}{2h} \sum_{i=1,3,\ldots}^{2N-1}

This formula makes it easy to successively increase the resolution until the required accuracy is achieved.

Finally I'd like to comment on the behaviour of the error when N gets really large. In the figure above you can see that the error levels off at about 2 × 10 − 7 when N increases above 2000. This error is above the numerical accuracy, which is about 10 − 8 in my calculations. When N increases further the error starts to increase again until, in the end it is even larger than before. The reason for this large error is that the individual terms in the sum become very small. Adding many small values up results in a loss of accuracy due to rounding errors.

1 comment:

  1. Nicely worked with problem, you are good in teaching about trapezium and you made your blog more attractive by solving each problem in easy way.
    free online math tutor chat