The Maths Geek: numerics

Showing posts with label numerics. Show all posts

Friday, 27 May 2011

Tutorial: Midpoint rule

In a previous post I spoke about the trapezium rule for integrating functions numerically. The trapezium rule could find the value of an integral with second order accuracy. In this post I will present an alternative method, the midpoint rule. The geometric interpretation of the midpoint rule appears to be more crude than the trapezium rule. Nonetheless, it turns out that the midpoint rule also has second order accuracy.

Like the trapezium rule, we divide our interval $[a; b]$ into $N$ equal sections of width $h = (b - a) / N$ . The area of each section is now approximated by the area of the rectangle of width h and whose height is given by the function value at the centre of the section. The centre of the nth section is given by
$x n + 1 / 2 = a + (n + 1 / 2) h$ .

Here $n = 0,1,..., N - 1$ . As before, we define $f n + 1 / 2 = f (x n + 1 / 2)$ . The area of the $n$ th rectangle is then given by $h f n + 1 / 2$ . And the value of the integral is given by

$\int_a^b f(x)\;dx = h( f_{1/2} + f_{3/2} + \ldots + f_{N-1/2}) = h\sum_{k=0}^{N-1}f_{k + 1/2}$

This is the final form for the midpoint rule. As with the trapezium rule, the accuracy of the integral can be increased by increasing N, i.e. decreasing the step size h. Again we can plot the error of the numerical integral against the number of steps on a logarithmic scale.

In figure the error of the integral

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

when evaluated with the midpoint and the trapezium rule is plotted. Clearly, the midpoint rule has a second order error, as has the trapezium rule. At first sight this is somewhat surprising since the midpoint rule makes no attempt to incorporate the change in the function value over the sections into the approximation.

To understand the accuracy of the approximation, we have to re-interpret the midpoint formula

h f n + 1 / 2

as the area of a trapezoid that intersects the function graph at the midpoint an which has a slope identical to the derivative of the function at that point. It is a simple, but very useful, geometric fact that the area of a trapezoid is the same as the area of a rectangle whose height is the average of the trapezoid's height.

Looking at the error curves a bit more closely, we can observe that the error of the midpoint rule is always slightly smaller than that of the trapezium rule. In fact when we divide the two error curves in the region where they decrease with N^2, we get an almost constant error ratio of about 1/2. This will lead us to improve the two methods. But this will be the subject of another post.

Wednesday, 11 May 2011

Successive Interval Halving

In a previous post I talked about the trapezium rule for approximating integrals. One nice feature was that one can successively increase the resolution without having to evaluate the function at points at which it has already been evaluated once. The formula I gave was

If one uses this method from the onset, this will result at the function being evaluated at the following points

Step $i$	$y = (x - a) / (b - a)$
0	0
1	1
2	1/2
3	1/4
4	3/4
5	1/8
6	3/8
7	5/8
8	7/8

and so on. A few days ago I was creating a demonstration and I wanted to create a longer table with these values. I didn't want to create it by hand so I was looking for a closed formula that produces the values in the second column from those in the first column. After some pondering over the problem I came up with the following formula

$y=(2n+1) 2^{-\lfloor\log_2(2n+1)\rfloor} - 1$

Here the $\lfloor.\rfloor$ denotes the floor function, i.e. the largest integer that is smaller or equal the argument. This formula will produce all the numbers in the table, including the zero, but not the one. The one, therefore, still has to be added by hand. Usually I would advise against using logarithms and exponentials in numerical algorithms but, because they are to the base 2, they can be implemented by bit shifting operations.

There is actually an interesting relation between $2 n + 1$ and $y$ when we write them as binary numbers. In fact, if we write $2 n + 1$ as binary and place a decimal point after the leading 1 we get $y + 1$ . Here is an example

2n+1	$100111012$
y+1	$1.0011101 2$
y	$0.0011101 2$

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

x i

between the integration boundaries

a

and

b

x i = a + h i

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

T N

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 = 1 e 6

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

2 N

, one does not have to evaluate the function at all the

2 N

points. Instead one can use the result

T N

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 \times 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.