Some basic aspects of numerical computation, including errors, interpolation, polynomials, integration, and some of the pitfalls: badly-conditioned problems and poor scaling. We use Maple and MATLAB to demonstrate these techniques.

High Performance Computing (HPC) is an engineering discipline which develops the technology (hardware, algorithms, and software) to deliver cost-effective computational results.

Allows you to choose a solution technique for your problem which is

• efficient
• stable

The question "why bother?" often arises, to which the answer is "to avoid investing time and effort in solving the wrong problem."

References names used in the text are shown in square brackets after the book name e.g. [NR] for Numerical Recipes.

• Abramowitz MA and Stegun IA (1971) Handbook of Mathematical Functions. (New York: Dover).
• Acton, F.S. (1990, corrected edition) Numerical Methods that work (Washington: Mathematical Association of America).
• Gander, W. and Hrebicek, J. (1991) Solving Problems in Scientic Computing using Maple and Matlab, (New York: Springer-Verlag). [GH]
• Gradshteyn, I.S. and Ryzhik, I.M. (1965, check latest edition) Table of Integrals, Series and Products (New York: Academic Press).
• Stoer, J. and Bulirsch R. (1993, second edition) Introduction to Numerical Analysis, (New York: Springer-Verlag).
• Press, W.H., Teulkolsky S.A., Vetterling, W.T. and Flannery B.P. (1992, second edition) Numerical Recipes in FORTRAN 77, (Cambridge: Cambridge University Press). See also the Fortran 90 and C versions. [NR]

There are several kinds of error

As with decimal (e.g. 3.) it is not possible to represent every number exactly in binary due to the use of a finite number of bits.

A floating point number is represented in the form

In MATLAB try the following:

EDU» a = 1.000000001^10

a =

1.00000001000000

EDU» b=a^(1/10)

b =

1.00000000100000

EDU» a-b

ans =

9.000000744663339e-009

Raising a number to a power, n, and then taking the nth root may not return the same number!

Similarly

EDU» tan(tan(tan(atan(atan(atan(1e50))))))

ans =

1.633177872838384e+016

In just three calculations, we were > 30 orders of magnitude out. Just because an equation can be algebraically simplified, it may not cancel numerically.

In the first case the computer was unable to represent our 10th root sufficiently accurately. In the second case, we chose a function with a singularity, so any small rounding errors in the computation would be magnified.

• Underflow and Overflow: a number may not be representable within the mantissa range allowed by the computer/ language. Try rescaling variables during the calculation.
• Machine addition is not associative:

In general this is because the lowest significant digits are lost. In general this can be addressed by adding numbers starting with the smallest first, or adding in pairs.

• Machine subtraction of two nearly equal numbers can cause a loss of accuracy
• Multiplication. The product of two m digit numbers is at most a 2m digit number: some machines represent multiplications using 2m digits to avoid loss of accuracy.
• Dividing a large number by a small number can cause overflow, and (/or) a loss of accuracy.

These rely on discretizing the integral and approximating the area under the curve by finite sums. In effect we replace the function by a suitable interpolating polynomial and integrate this.

Open formulas do not use the value of the function at its endpoints (a and b). An open formula can be useful if the function has an integrable singularity at one end.

Closed formula use the value of the functions at its endpoints. The most widely known of these is the trapezium rule, which is equivalent to performing linear interpolation between the points.

EDU» fplot('humps',[-1,2],'r')

EDU» x = linspace(-1, 2, 30);

EDU» y = humps(x);

EDU» hold on

EDU» stem(x,y,':')

EDU» plot(x,y,':')

EDU» title('Trapezium Rule for Humps (30 points)')

EDU» hold off

MATLAB allows the area to be found using:

trapz(x,y)

ans =

26.2102

The error in the trapezium rule depends upon the value of the derivatives at the end points, so if these derivatives vanish, or if the function is being integrated over one period, the trapezium rule gives good value for money.

From the plot, it is clear that some of the area has been missed. We could add more points along the range of interest, until the approximation was good enough, however this would result in a considerable amount of extra work. There are two solutions to this:

In the previous cases, we picked equally spaced abscissas. In the case of the peak in our function, simple linear interpolation missed some of the area. Why not exploit our knowledge of the function to pick unequal spacing for the abscissa to better reflect where the function has most of its area? We must then weight each of our abscissa correctly

The derivation of the formula and weights relies on the expansion of a function in terms of suitable orthogonal polynomials. However, there are many sources of information on Gaussian Integration weights and abscissas (e.g. section 4.5 of NR for more details and Abramowitz and Stegun Chapter 25). Naturally, you can use Maple to derive weights and abscissa not listed in tables- efficient procedures exist to do this (see [GH] for details).

These are harder, due to the large number of samples of the function required in order to achieve a given accuracy.

• Use any available symmetry of the system to reduce the region of integration to a minimum.
• For N-dimensional integration where N is large of an arbitrary function, there are good Monte Carlo methods. Some of these are very sophisticated and focus on regions of interest (i.e. where the integral is large) to reduce the errors of this approximation.

Evaluation of a polynomial should be performed by a using a Horner scheme. Consider the following example in Maple:

> y := a*x^4 + b*x^3 + c*x^2 + d*x + e;

> convert(y,'horner',x);

>

e + (d + (c + (b + a x) x) x) x

Why? The nested Horner representation takes only O(N) multiplications, whereas the conventional notation requires O(N²) multiplications. By storing the coefficients of the polynomial in an array, the Horner scheme allows a succinct way to evaluate the polynomial. This method also allows one to evaluate a polynomial and its derivatives in a single loop (such as one might require for a Newton-Raphson iteration to find a root)

There are clever ways to:

• Multiply two polynomials. The product of two polynomials can be calculated by taking the Fourier transform of the polynomial coefficients, multiplying, and then taking the inverse Fourier transform (an example of convolution).
• Divide a polynomial of degree (n - 1) by a factor (x - a), without performing any division. This is "synthetic division" and can be extended to divide to general polynomials. If you find one root of an order N polynomial, it can be used to divide out this root and work with the remaining polynomial of degree N-1 (this process is called "deflation", and is also used in eigenvalue problems), but beware of rounding errors.

There is considerable literature about finding the roots of polynomials. The problems involved with such an innocuous task are useful to motivate the study of numerical analysis. Problems include:

• A general polynomial may have complex zeros.
• Conventional methods may be fooled by equations with multiple roots.
• Newton's method may diverge
• Finding roots of a polynomial is ill-conditioned. A small change in the coefficients can dramatically affect the solution. Consider the following example using Maple

>solve(x^4+4*x^3+6*x^2+4*x+1,x);

-1, -1, -1, -1

Now change the coefficient in front of the x term:

>solve(x^4+4*x^3+6*x^2+4.01*x+1,x);

-1.340248595,

-.9750019523 - .3172265943 I,

-.9750019523 + .3172265943 I,

-.7097475007

and MATLAB:

EDU» roots([1,4,6,4.01,1])

ans =

-1.3402

-0.9750 + 0.3172i

-0.9750 - 0.3172i

-0.7097

So beware. In fact this is example is pathologically ill-conditioned, since the original function has a quadruple root at x = -1.

EDU» A = gallery(3)

A =

-149 -50 -154

537 180 546

EDU» rh1 = [-711; 2535; -120], x = A \ rh1

rh1 =

-711

2535

-120

x =

1.0000

2.0000

3.0000

EDU» rh2 = rh1 + [1; 1; 0.1], x = A \ rh2

rh2 =

-710

2536

-120

x =

99.6667

-312.0667

9.5000

What has gone wrong? We perturbed the right hand side by 0.1%, well within the bounds of experimental error, but the solution changed by two orders of magnitude in one case! The matrix is very badly conditioned. We can measure this using:

EDU» cond(A)

ans =

A different problem that can arise is "poor scaling", in which matrix elements vary over several orders of magnitude. Numerical rounding errors cause a loss of accuracy in the final solution. Fortunately MATLAB's built-in functions use techniques such as pivoting to reduce the effects of poor scaling. If you find a difference between MATLAB's results and those from a numerical FORTRAN or C routine, you should check to see whether the coefficients of the problem are poorly scaled. MATLAB will have attempted to control any loss of accuracy.

Before applying a numerical method to a problem understand carefully any possible problems which may arise, and ensure that the implementation you have chosen is

• efficient
• stable

In the next numerical analysis lectures we will consider the problems involved in more advanced numerical methods, such as the solution of ordinary and partial differential equations, finding eigenvalues, and optimization methods (such as genetic algorithms). A number of specific numerical problems arise in parallel programming, including use of halos in computation, parallel linear algebra, long range forces, and random number generation.