152
6 Computing Integrals and Testing Code
criterion is that the exact solution is known. The other criterion is that the numerical
algorithm should produce the exact solution, within machine precision, for any
chosen n.
The second criterion may seem strange at first, but take the trapezoidal method,
for example, and consider an integral like
b
a (6x − 4)dx. The integrand is here a
straight line, and if you sketch it in a coordinate system along with some relevant
trapezoids, you realize that the topmost side of each trapezoid comes exactly on the
straight line. This is the case for whatever number n of trapezoids you may choose to
use. Thus, the trapezoidal method should solve that problem without any numerical
error. 4 We can therefore pick some linear function and construct a test function that
checks for equality between the exact solution and the numerical approximation
produced by our implementation of the trapezoidal method.
Note that, when testing, numbers should not just be taken out of the air. For
example, a specific test case can be
4.4
1.2 (6x − 4)dx. This integral involves an
“arbitrary” interval [1.2, 4.4] and an “arbitrary” linear function f (x) = 6x − 4. By
“arbitrary”, we mean expressions where the special numbers 0 and 1 are avoided,
since these have special properties in arithmetic operations (e.g., forgetting to
multiply is equivalent to multiplying by 1, and forgetting to add is equivalent to
adding 0).
Demonstrating Correct Convergence Rates Also for these unit tests we choose
problems for which the exact solution is known. However, contrary to the previous
test procedure, we now work with problems for which the numerical algorithm does
not give zero approximation error. Normally, unit tests must be based on that kind
of problems. Thus, the answer we get from our code will contain an approximation
error, and since we know the exact solution, we may compute the size of this error.
Unfortunately, this is not too helpful, since we have little chance telling if this error
is what we should have got for the particular n used!
Now the convergence rate comes in handy. If we know (or have reason to assume)
that the numerical error has a certain asymptotic behavior when n → ∞, we know
at what rate the numerical error should be reduced. The idea of a corresponding unit
test is then to run the algorithm for some n values, compute the error (the absolute
value of the difference between the exact solution and the one produced by the
numerical method), and check that the error has approximately correct asymptotic
behavior. For the trapezoidal and midpoint methods in particular, this means that the
error should become proportional to n −2 when n → ∞.
Let us develop a more precise method for such unit tests based on convergence
rates. Consider a set of q + 1 experiments with various n: n 0 , n 1 , n 2 , . . . , n q . We
compute the corresponding errors E 0 , . . . , E q . For two consecutive experiments,
number i and i − 1, we have the error model
E i = Cn
−r
i ,
(6.26)
E i−1 = Cn
−r
i−1 .
(6.27)
4 In fact, so would the midpoint method! This is because, for each rectangle, the error to each side
of the midpoint is equally large with opposite signs, meaning that they cancel each other.
6 Computing Integrals and Testing Code
criterion is that the exact solution is known. The other criterion is that the numerical
algorithm should produce the exact solution, within machine precision, for any
chosen n.
The second criterion may seem strange at first, but take the trapezoidal method,
for example, and consider an integral like
b
a (6x − 4)dx. The integrand is here a
straight line, and if you sketch it in a coordinate system along with some relevant
trapezoids, you realize that the topmost side of each trapezoid comes exactly on the
straight line. This is the case for whatever number n of trapezoids you may choose to
use. Thus, the trapezoidal method should solve that problem without any numerical
error. 4 We can therefore pick some linear function and construct a test function that
checks for equality between the exact solution and the numerical approximation
produced by our implementation of the trapezoidal method.
Note that, when testing, numbers should not just be taken out of the air. For
example, a specific test case can be
4.4
1.2 (6x − 4)dx. This integral involves an
“arbitrary” interval [1.2, 4.4] and an “arbitrary” linear function f (x) = 6x − 4. By
“arbitrary”, we mean expressions where the special numbers 0 and 1 are avoided,
since these have special properties in arithmetic operations (e.g., forgetting to
multiply is equivalent to multiplying by 1, and forgetting to add is equivalent to
adding 0).
Demonstrating Correct Convergence Rates Also for these unit tests we choose
problems for which the exact solution is known. However, contrary to the previous
test procedure, we now work with problems for which the numerical algorithm does
not give zero approximation error. Normally, unit tests must be based on that kind
of problems. Thus, the answer we get from our code will contain an approximation
error, and since we know the exact solution, we may compute the size of this error.
Unfortunately, this is not too helpful, since we have little chance telling if this error
is what we should have got for the particular n used!
Now the convergence rate comes in handy. If we know (or have reason to assume)
that the numerical error has a certain asymptotic behavior when n → ∞, we know
at what rate the numerical error should be reduced. The idea of a corresponding unit
test is then to run the algorithm for some n values, compute the error (the absolute
value of the difference between the exact solution and the one produced by the
numerical method), and check that the error has approximately correct asymptotic
behavior. For the trapezoidal and midpoint methods in particular, this means that the
error should become proportional to n −2 when n → ∞.
Let us develop a more precise method for such unit tests based on convergence
rates. Consider a set of q + 1 experiments with various n: n 0 , n 1 , n 2 , . . . , n q . We
compute the corresponding errors E 0 , . . . , E q . For two consecutive experiments,
number i and i − 1, we have the error model
E i = Cn
−r
i ,
(6.26)
E i−1 = Cn
−r
i−1 .
(6.27)
4 In fact, so would the midpoint method! This is because, for each rectangle, the error to each side
of the midpoint is equally large with opposite signs, meaning that they cancel each other.