The series of goals above are briefly stated, but illustrated here in detail for a special, simple case study: numerical integration by the Trapezoidal rule.
Many science courses now have examples and exercises involving implementation and application of numerical methods. How to structure and verify such numerical programs has, unfortunately, received little attention in university education and the literature. Students and teachers occasionally write programs that are too tailored to the problem at hand instead of being a good starting point for future extensions, and testing is often limited to running a case where the answer seems reasonable. The standards of computing need to be raised to the levels found in experimental physics, chemistry, and biology.
Integrate the function \( g(t)=\exp{(-t^4)} \) from -2 to 2 using the Trapezoidal rule, defined by $$ \begin{equation} \int_a^b f(x)dx \approx h\left( {1\over2}(f(a) + f(b)) + \sum_{i=1}^{n-1} f(a+ih)\right), \quad h = (b-a)/n \tag{1} \end{equation} $$
Many will attempt to solve the problem by this simple program in Matlab:
a = -2; b = 2;
n = 1000;
h = (b-a)/n;
s = 0.5*(exp(-a^4) + exp(-b^4));
for i = 1:n-1
s = s + exp(-(a+i*h)^4);
end
r = h*s;
r
The solution is minimalistic and correct. Nevertheless, this solution has a common pedagogical and software engineering flaw: a special function \( \exp(-t^4) \) is merged into a general algorithm (1) for integrating an arbitrary function \( f(x) \).
The writer of the program runs it and reports the result: 1.81280494737. How can one assess that this result is correct? There is no exactly known result to compare with. Also, the program above is not well suited for switching to an integrand where we can compare with an exact answer, because several lines need modification.
A fundamental software engineering practice is to use functions for splitting a program into natural pieces, and if possible, make these functions sufficiently general to be reused in other problems. In the present problem we should strive for the following principles:
Trapezoidal.m containing
function r = Trapezoidal(f, a, b, n)
% Numerical integration from a to b
% with n intervals by the Trapezoidal rule
f = fcnchk(f);
h = (b-a)/n;
s = 0.5*(f(a) + f(b));
for i = 1:n-1
s = s + f(a+i*h);
end
r = h*s;
The special \( g(t) \) function can be implemented in a separate file g.m
or put in the main program. The function becomes
function v = g(t)
v = exp(-t^4)
end
Finally, a specialized main program (main.m) solves the problem at hand:
a = -2; b = 2;
n = 1000;
result = Trapezoidal(@g, a, b, n);
disp(result);
exit
The important feature of this solution is that Trapezoidal.m can be
reused for "any" integral. In particular, it is straightforward
to also integrate an integrand where we know the exact result.
An advantage of having the \( g(t) \) as a separate function is that we can easily send this function to a different integration method, e.g., Simpson's rule.
Both Solution 1 and Solution 2 are readily implemented in Python.
However, functions in Python do not need to be located in separate
files to be reusable, and therefore there is no psychological barrier
to put a piece of code inside a function. The consequence is that a
Python programmer is more likely to go for Solution 2.
(This may be the reason why the author has observed
scientific Python codes to have better design
than Matlab codes - modularization comes more natural.)
The relevant
code can be placed in a single file, say main.py, looking as
follows:
def Trapezoidal(f, a, b, n):
h = (b-a)/float(n)
s = 0.5*(f(a) + f(b))
for i in range(0,n,1):
s = s + f(a + i*h)
return h*s
from math import exp # or from math import *
def g(t):
return exp(-t**4)
a = -2; b = 2
n = 1000
result = Trapezoidal(g, a, b, n)
print result
This solution acknowledges the fact that the implementation is a generally applicable function, just as the Trapezoidal formula.
However, a small adjustment of this implementation will make it
much better. If somebody wants to reuse the Trapezoidal function
for another integral, they can import Trapezoidal from the
main.py file, but the special problem will be executed as part of
the import. This is not desired behavior when solving another problem.
Instead, our special exercise problem should be isolated in
its own function and called from a test block in the file (to avoid
being executed as part of an import).
This is the general software design of modules in Python.
We therefore rewrite the code in a new file Trapezoidal.py:
def Trapezoidal(f, a, b, n):
h = (b-a)/float(n)
s = 0.5*(f(a) + f(b))
for i in range(1,n,1):
s = s + f(a + i*h)
return h*s
def _my_special_problem():
from math import exp
def g(t):
return exp(-t**4)
a = -2; b = 2
n = 1000
result = Trapezoidal(g, a, b, n)
print result
if __name__ == '__main__':
_my_special_problem()
Now we have obtained the following important features:
Trapezoidal.py is a module offering the widely applicable
function Trapezoidal for integrating "any" function.Trapezoidal.py is run as a program, the if test is true and
the special integral of \( g \) is computed.from Trapezoidal import Trapezoidal, the if test
is false and nothing gets computed.An integral part of any implementation is verification, i.e., to bring evidence that the program works correctly from a mathematical point of view. (A related term, validation, refers to bringing evidence that the program produces results in accordance with observations of nature, but this is not of direct interest in this integration context.)
The intuitive approach to testing is to compare results of a program with known mathematical results. For example, we can choose some function, say \( \sin t \), and differentiate it to obtain an integrand that we can easily integrate by hand and thereby get a precise number for the integral. Integrating \( \int_{-2}^2 \cos t\,dt \) gives the exact result 1.81859485365. The program with the Trapezoidal rule reports 1.81859242886, so the error \( 2.42\cdot 10^{-6} \) is "small". However, we have no idea if this error is just the approximation error in the numerical method or if the program has a bug too! What if the error were \( 1.67\cdot 10^{-3} \)? It is impossible to say whether this answer is the correct numerical result or not. Actually, this error contains both the approximation error and a bug where the loop goes over \( 0,1,\ldots,n-1 \).
So, comparison of a numerical approximation with an exact answer does not say much unless the error is "huge" and therefore clearly points to fundamental bugs in the code.
For most numerical methods there are only two good verification methods:
The Trapezoidal rule is obviously exact for linear integrands. Therefore, we should test an "arbitrary" linear function and check that the error is close to machine precision. This is done in a separate function in a separate file test_Trapezoidal.py:
from Trapezoidal import Trapezoidal
from Trapezoidal_vec import Trapezoidal as Trapezoidal_vec
def linear():
"""Test linear integrand: exact result for any n."""
def f(x):
return 8*x + 6
def F(x):
"""Anti-derivative of f(x)."""
return 4*x**2 + 6*x
a = 2
b = 6
exact = F(b) - F(a)
numerical = Trapezoidal(f, a, b, n=4)
error = exact - numerical
print '%.16f' % error
The output of calling linear()
is in this case zero exactly, but in general one must expect
some small rounding errors in the numerical and exact result.
The function linear performs the test, but it would be better to
integrate the test into a testing framework such that we with
one command can execute a comprehensive set of tests. This makes it
easy to run all tests after every small change of the software.
Students should adopt such compulsory habits from the software industry.
The dominating type of test frameworks today is based on what is called unit testing in software engineering. It means that we pick a unit in the software and write a function (or class) that runs the test after certain specifications:
test_.AssertionError exception (in Python) is raised, otherwise the function runs silently.
In our case, a proper test function means the following rewrite
of the function linear:
def test_linear():
"""Test linear integrand: exact result for any n."""
def f(x):
return 8*x + 6
def F(x):
"""Anti-derivative of f(x)."""
return 4*x**2 + 6*x
a = 2
b = 6
expected = F(b) - F(a)
tol = 1E-14
computed = Trapezoidal(f, a, b, n=4)
error = abs(expected - computed)
msg = 'Trapezoidal: expected=%g, computed=%g, error=%g' % \
(expected, computed, error)
assert error < tol, msg
The code is basically the same, but we comply to the rules in the three
bullet points above. The assert
statement has the test as error < tol, with msg as an optional
message that is printed only if the test fails (error < tol is False).
The msg string can be left out and it suffices to do assert error < tol.
The reason why we comply to testing frameworks is that we can use software
like nose or pytest to automatically find all our tests and execute them.
We put tests in files or directories starting with test and run
one of the commands
Terminal> nosetests -s -v .
Terminal> py.test -s -v .
All functions with names test_*() in all files test*.py in all
subdirectories with names test* will
be run, and statistics about how many tests that failed will be printed.
The tests should be run after every modification of the software.
We integrated by hand the linear function used in the test above. In more complicated cases it would be safer to use symbolic computing software to carry out the mathematics. Here we demonstrate how to use the Python package SymPy to do the integration:
def test_linear_symbolic():
"""Test linear integrand: exact result for any n."""
import sympy as sym
# Define a linear expression and integrate it
x = sym.symbols('x')
f = 8*x + 6 # Integrand for this test
F = sym.integrate(f, x)
# Verify symbolic computation: F'(x) == f(x)
assert sym.diff(F, x) == f
# Transform expressions f and F to Python functions of x
f = sym.lambdify([x], f, modules='numpy')
F = sym.lambdify([x], F, modules='numpy')
# Run one test with fixed a, b, n, for scalar and
# vectorized implementation
a = 2
b = 6
expected = F(b) - F(a)
tol = 1E-14
for func in Trapezoidal, Trapezoidal_vec:
computed = func(f, a, b, n=4)
error = abs(expected - computed)
msg = 'expected=%g, computed=%g, error=%g' % \
(expected, computed, error)
assert error < tol, msg
Note that we now also test both the scalar and the vectorized implementations
of the Trapezoidal rule (see the section High-performance computing: vectorization for explanation of the
vectorized version Trapezoidal_vec).
It is easy in Python to loop over functions (with a variable like func).
We could also just compare the result of Trapezoidal_vec to that of
Trapezoidal when the latter is verified against the expected value.
Let us change the integration limits in our test example to \( a=2\cdot 10^8 \) and \( b=6\cdot 10^9 \). The computed error in this case is 16384 (!). Hence the tolerance must be set to (e.g.) \( 2\cdot 10^5 \) (!). In general, the tolerance depends on the magnitude of the numbers involved in the computations. To avoid this dependence, one should use relative errors:
error = abs(expected - computed)/abs(expected)
Now, a tolerance of \( 10^{-14} \) is adequate for the test even if
the numbers expected and computed are large.
A function test_linear_symbolic_large_limits in the
file test_Trapezoidal.py is a test function for a case with large limits
and use of the relative error.
Let us extend the verification with a case where we know the exact answer of the integral, but we do not know the approximation error. The only knowledge we usually have about the approximation error is of asymptotic type. For example, for the Trapezoidal rule we have an expression for the error from numerical analysis: $$ E = - \frac{(b-a)^3}{12n^2}f''(\xi),\quad\xi\in [a,b]. $$ Since we do not know \( \xi \), which is some number in \( [a,b] \), we cannot compute \( E \). However, we realize that the error has an asymptotic behavior as \( n^{-2} \): $$ E = Cn^{-2}, $$ for some unknown constant \( C \). If we compute two or more errors for different \( n \) values, we can check that the error decays as \( n^{-2} \). In general, when verifying the implementation of a numerical method with discretization parameter \( n \), we write \( E=Cn^r \), estimate \( r \), and compare with the exact result (here \( n=-2 \)).
More precisely, we perform a set of experiments for \( n=n_0, n_1, \ldots, n_m \), where we empirically estimate \( r \) from two consecutive experiments: $$ \begin{align*} E_i &= Cn_i^r,\\ E_{i+1} &= Cn_{i+1}^r. \end{align*} $$ Dividing the equations and solving with respect to \( r \) gives $$ r = \frac{\ln(E_i/E_{i+1})}{\ln(r_i/r_{i+1})}$$ As \( i=0,\ldots,m-1 \), the \( r \) values should approach the value \( -2 \).
It is easy to use the method of manufactured solutions to construct a test problem. That is, we first choose the solution, say the integral is given by \( F(b)-F(a) \), where $$ F(x) = e^{-x}\sin (2x).$$ Then we fit the problem to accept this solution. In the present case it means that the integrand must be \( f(x)=F'(x) \). We use for safety symbolic software to calculate \( f(x) \). Thereafter, we run a series of experiments where \( n \) is varied, we compute the corresponding convergence rates \( r \) from two consecutive experiments and test if the final \( r \), corresponding to the two largest \( n \) values, is sufficiently close to the expected convergence rate \( -2 \):
def test_convergence_rate():
import sympy as sym
# Construct test problem
x = sym.symbols('x')
F = sym.exp(-x)*sym.sin(2*x) # Anti-derivative
f = sym.diff(F, x) # Integrand for this test
# Turn to Python functions
f = sym.lambdify([x], f, modules='numpy')
F = sym.lambdify([x], F, modules='numpy')
a = 0.1
b = 0.9
expected = F(b) - F(a)
# Run experiments (double n in each experiment)
n = 1
errors = []
for k in range(28):
n *= 2
computed = Trapezoidal(f, a, b, n)
error = abs(expected - computed)
errors.append((n, error))
print k, n, error
# Compute empirical convergence rates
from math import log as ln
estimator = lambda E1, E2, n1, n2: ln(E1/E2)/ln(float(n1)/n2)
r = []
for i in range(len(errors)-1):
n1, E1 = errors[i]
n2, E2 = errors[i+1]
r.append(estimator(E1, E2, n1, n2))
expected = -2
computed = r[-1] # The "most" asymptotic value
error = abs(expected - computed)
tol = 1E-3
msg = 'Convergence rates: %s' % r
print errors
assert error < tol, msg
The empirical convergence rates are in this example
-2.022, -2.0056, -2.0014, -2.00035, -2.000086, -2.000022,
-2.0000054, -2.0000013, -2.00000033
Although the rates are known to approach \( -2 \) as \( n\rightarrow\infty \), the rates are close to \( -2 \) even for large \( n \) (such as \( n=4 \)). A rough tolerance is often used for convergence rates, for instance 0.1, but here we may use a smaller one if desired.
Known analytical solutions are of value in convergence rate tests, but if they are not available, or restricted to very simplified cases, the method of manufactured solutions, where we solve a perturbed problem fitted to a constructed exact solution, is also a very useful technique.
In Matlab, one must decide whether to use a class-based system for unit
testing or just write test functions that mimic the behavior of the
Python test functions for the nose and pytest frameworks. Here is an
example on doing the test_linear() function in Matlab:
function test_trapezoidal_linear
%% Check that linear functions are integrated exactly
f = @(x) 8*x + 6;
F = @(x) 4*x**2 + 6*x; %% Anti-derivative
a = 2;
b = 6;
expected = F(b) - F(a);
tol = 1E-14;
computed = trapezoidal(f, a, b, 4);
error = abs(expected - computed);
assert(error < tol, 'n=%d, error=%g', n, error);
end
end
test_trapezoidal_linear()
There is, unfortunately, no software available to run all tests in all files in all subdirectories and report on the success/failure statistics, but it is quite straightforward to write such software.
Verification in terms of measuring convergence rates usually gives a very good insight into approximation errors, but the verification results may be affected by rounding errors, depending on the type of algorithm. For the scalar implementation of the Trapezoidal rule, rounding errors start to affect the results around \( n=2^{24}\approx 16 \) million points. Other algorithms are much more sensitive to rounding errors. For example, a numerical derivative like $$ f''(x)\approx \frac{f(x+h)-2f(x) + f(x+h)}{h^2},$$ may be subject to rounding errors for moderate values of \( h \). Here is an example with \( f(x)=x^{-6} \). An exact answer is \( f''(1)=42 \), but numerical experiments for with \( h=10^{-k} \) for various \( k \) values end up with
| \( k \) | numerical \( f'' \) |
| 1 | 44.61504 |
| 2 | 42.02521 |
| 3 | 42.00025 |
| 4 | 42.00000 |
| 5 | 41.99999 |
| 6 | 42.00074 |
| 7 | 41.94423 |
| 8 | 47.73959 |
| 9 | -666.13381 |
| 10 | 0.00000 |
| 11 | 0.00000 |
| 12 | -666133814.8 |
| 13 | 66613381477.5 |
| 14 | 0.00000 |
The error starts to increase rather than decrease for \( h>10^{-5} \), and this is because the rounding error is (much) bigger than the approximation error in the formula.
We discuss here how some of the learning outcomes from the section General learning outcomes for computing competence can be incorporated in the exercise with the Trapezoidal rule. We restrict programming examples to use Python.
This author has seen a lot of programs used for teaching which apply vectorization without explicit notice. Vectorization is a technique in high-level languages like IDL, MATLAB, Python, and R for removing loops and speed up computations. Unfortunately, the "distance" from the mathematical algorithm to vectorized code is larger than to a plain loop as we used above. Vectorization therefore tends to confuse students who are not well educated in the techniques. For example, the Trapezoidal rule can be vectorized as
import numpy as np
def Trapezoidal(f, a, b, n):
x = np.linspace(a, b, n+1)
return (b-a)/float(n)*(np.sum(f(x)) - 0.5*(f(a) + f(b)))
The code is correct, but it takes some thinking to realize why these
lines compute the formula (1). Because of the sum function,
we need to adjust the summation result such that
the weight of the end points becomes correct.
The scalar implementation of the Trapezoidal rule computes one f(x)
at the time and uses very little memory, actually only 4 float
variables. The vectorized version, however, computes the
function values at all points x (\( n+1 \) float objects)
at once and therefore requires the storage of about \( n \) float objects.
This is a significant difference between the vectorized and scalar
versions. The vectorized version may run out of memory if we
want very accurate results and hence a large \( n \).
An important observation for parallelization of the Trapezoidal
rule is that all the function evaluations are independent of each
other so these can be performed in parallel. Typically, with \( m \)
compute units we can distribute int(n/m) function evaluations
to the first \( m-1 \) units and int(n/m) + n % m to the last one.
Each unit must compute the sum of the evaluations and
communicate to one master unit or to all other units.
The master or all units must then sum up all the partial sums,
subtract \( \frac{1}{2}(f(a) + f(b)) \) and multiply by \( h \) to get the
final answer.
Vectorized algorithms often lend themselves to automatic parallelization. In fact, the Numba tool can automatically parallelize Numerical Python code. Looking at the vectorized Trapezoidal implementation
x = np.linspace(a, b, n+1)
I = (b-a)/float(n)*(np.sum(f(x)) - 0.5*(f(a) + f(b)))
and assuming that n is large, we realize that np.linspace must
create the vector on \( m \) compute units, each with its own memory.
In the next expression, f(x) leads to application of f on the
piece of x that is on each compute unit. Then np.sum creates
partial sums of f(x) on each compute unit and distributes the
results to all other units. No more (distributed) vectors are involved,
so the remaining scalar operations can be carried out on every unit,
and the final result of the integral is then available on each individual
compute unit.
We think it is fundamental that such reasoning is well known among students. Traditionally, thinking about parallelism has not been in focus unless also demanding technical implementations in terms of MPI is also taught. However, laptops will soon be powerful parallel computing platforms, so knowing how to write code that lend itself to easy parallelization by tools such as NumPy and Numba is key. How parallel code is actually implemented may be pushed to a more specialized courses.
Since the asymptotic behavior of approximation errors is so fundamental for the most common verification technique (i.e., checking convergence rates), students should be well motivated for diving more into the mathematics behind the various formulas they use in test functions.
The Trapezoidal rule is primarily a pedagogical tool for obtaining a good understanding numerical integration and what integration is. For professional use, one will apply more sophisticated algorithms, for instance algorithms that deliver an estimate of the integral with a specified error tolerance.
In the scientific Python eco system we have the quad method
for sophisticated integration, from
the famous QUADPACK Fortran library and made available in the SciPy package.
Here we integrate \( \int_{-2}^2\cos t\,dt \) with a relative error
tolerance of \( 10^{-12} \):
>>> import scipy.integrate
>>> from math import cos
>>> I, error = scipy.integrate.quad(cos, -2, 2, epsrel=1E-12)
>>> I
1.8185948536513632
>>> error
2.4124935390890847e-14
We want to compute \( I=\int_a^{2}\cos t\,dt \), but \( a \) is an uncertain parameter. Suppose \( a \) can be modeled as a normally distributed stochastic variable with mean -2 and standard deviation 0.2. What is the corresponding uncertainty in \( I \)? The simplest statistics reflecting the uncertainty of \( I \) is the mean and the standard deviation. Monte Carlo simulation is the simplest method for computing the uncertainty.
from Trapezoidal_vec import Trapezoidal
import numpy as np
N = 100000 # Monte Carlo samples
a = np.random.normal(loc=-2, scale=0.2, size=N) # N samples of a
I = np.zeros(N) # Responses (integrals)
for i in range(N):
I[i] = Trapezoidal(np.cos, a[i], 2, n=1000)
print 'Integral of cos(t) from t=-2 to t=2:', np.sin(2) - np.sin(-2)
print 'Mean value of uncertain integral:', np.mean(I)
print 'Standard deviation of uncertain integral:', np.std(I)
The output becomes
Integral of cos(t) from t=-2 to t=2: 1.81859485365
Mean value of uncertain integral: 1.80044133926
Standard deviation of uncertain integral: 0.0856614641125