Numerical Analysis

Numerical Analysis II

An overview of advanced numerical methods, such as solution of ordinary and partial differential equations, and optimization methods (such as genetic algorithms).

(To read this with no frames type this location into your web browser: training/numan_ii.html)

High Performance Computing

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

Numerical Analysis

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."

Bibliography

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

Acton, F.S. (1990, corrected edition) Numerical Methods that work (Washington: Mathematical Association of America).
Farlow, S.J. (1993) Partial Differential Equations for Scientists and Engineers (New York: Dover).
Gradshteyn, I.S. and Ryzhik, I.M. (1965, check latest edition) Table of Integrals, Series and Products (New York: Academic Press).
Hackbusch, W. (1985) Multi-Grid methods and applications, (New York: Springer-Verlag)
Iserles, A. (1996) A First Course in the Numerical Analysis of Differential Equations (Cambridge: Cambridge University Press). ISBN 0-521-556550-4
Michalewicz, Z. (1994, 2nd Edition) Genetic Algorithms + Data Structures = Evolution Programs (New York: Springer-Verlag). Note that there is now third edition of this book.
Nash, S.G., and Sofer, A. (1996) Linear and Nonlinear Programming (New York: McGraw-Hill)
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]
Stoer, J. and Bulirsch R. (1993, second edition) Introduction to Numerical Analysis, (New York: Springer-Verlag).
Zwillinger, D. (1997, third edition) Handbook of differential equations, (New York: Academic Press) ISBN 0-12-784396-5

See also Ken Thomas’ multigrid bibliography:

http://www.ecs.soton.ac.uk/~kst/mgrid.sem/bib.htm

Differential Equations

Most physical phenomena in science and engineering can be described in terms of an underlying differential equation (or set of differential equations).

The general strategy will always be to re-write differentials in terms of finite differences:

( 1 )

There are many different choices for finite difference schemes, and appropriate use of the correct scheme is at the heart of ensuring that the solution method is stable.

Ordinary Differential Equations

As with integration, a large number of simple differential equations can be solved analytically. If the answer yields a simple form, then it is best to work with this.

The general form is a set of coupled equations in terms of the first derivatives of a set of functions

e.g. y’ = f (t , y)

Higher order ODEs, with for example y’’, can be written in terms of a system of first order ODEs

Initial value problems

We know the y_i for some starting x_j, and we would like to know them for later x_j. Conceptually we start at the initial value and march out cautiously along the function. Methods to use are the Runge-Kutta technique (see Numerical Recipes 16.1), or adaptive techniques, which attempt to control the errors in the solution by taking smaller steps where necessary (see NR, 16.2).

Two point boundary values problems

Here the value of y is fixed at (usually) the beginning and the end of the region of interest. These are much harder than initial value problems, since we must only consider functions which obey the boundary conditions. Conceptually we start at the beginning of the region and "shoot out" to the end of the region: on successive iterations we aim to match the boundary conditions.

Stiffness

If the ratio of the largest : smallest eigenvalues of the system matrix is large, then the system of equations will be stiff. This ratio is measures by the condition number of the matrix. In this case, great care must be taken: solving stiff problems can lead to a dramatic loss of precision. See NR 16.6, or Stoer and Bullirsch 7.2.16. Consider using an implicit method (see below).

Implicit vs Explicit methods

In explicit methods, it is possible to solve (at a point) directly for the value of all unknowns in the finite difference scheme.

In an implicit scheme, there is no explicit formula at a points, only a set of simultaneous equations which must be solved to give the solution over the whole grid. Implicit methods are stable for all step sizes.

Partial Differential Equations

In contrast to ODEs, a partial differential equation relates a number of variables and their derivatives.

Hyperbolic Equations

Wave equation

( 2 )

Parabolic Equations

Diffusion Equation

( 3 )

Elliptic Equations

Poisson equation

( 4 )

Initial Value Problems

Here we wish to propagate the solution forward in time.

Boundary Value Problems

The solution is static.

Dirichlet conditions specify the boundary values around the region.
Von Neumann conditions specify the derivatives at the boundary. For Neumann conditions, some other information must be specified, or the solution will not be unique.

Stability

A reasonable looking scheme to solve a PDE may converge, but if the wrong finite difference scheme is used, then it may converge to the solution of the wrong problem! Stability is checked using, for example, Fourier analysis.

For diffusion and wave type equations, the ratio is between the time steps and space steps is constrained by the "Courant condition". See NR 19.1

The Multigrid method

This is a most important method for the solution of linear (and non-linear) elliptic equations. It is O(N), where N is the number of grid points. The key idea is to introduce a hierarchy of grids. Methods, such as Gauss-Seidel tend to stall, since they do not attenuate smooth error modes. However, on a coarser grid, it is not possible to represent these smooth errors, so by alternating between fine and coarse grids, the system relaxes more quickly.

Finite Element Methods

Used extensively in engineering for elliptic PDE solution over irregular bodies. There is extensive literature on these techniques and where they are useful. In principle a finite element calculation consists of

Discretizing the domain into simple shapes (usually triangles)
Assembling a stiffness and load matrix
Rewriting the differential equation as a variational principle
Constructing element stiffness matrices and load vectors element by element
Performing a minimization (usually of a sparse system)

Optimization Methods

The general optimization problem can be stated as follows:

"Minimize F(x) subject to g(x)"

F(x) is the objective function, and the set of functions g(x) are the constraints, which may be linear or non-linear.

Often the constraints are incorporated into the minimization by using a penalty function

Note that maximization of F(x) is the same as miminization of (-F(x))

This problem is one of the most widely solved on a day-to-day basis in science and engineering.

General Optimization Algorithm

Specify an initial guess, x₀, for the solution

for i = 0, 1, …

if x_i is optimal, then stop else form x_i+1, and repeat (3)

With derivatives or without?

In the simplest case where we wish to minimize F(x), the most important question is whether you can compute the derivatives of F(x).

With Derivatives

Conjugate gradient type methods, such as Fletcher-Reeves. These rely on a set of successive line minimizations to perform the optimization

Quasi Newton methods. These rely on computing an approximation to the underlying Hessian Matrix

No derivatives

The simplex method (Nelder-Mead simplex, not linear programming simplex). This a robust method to "get something working", and should not be overlooked.

Powell’s method, which derives information about the local gradient of the function, and is quadratically convergent.

If the function is very complicated, consider a probabilistic method such as simulated annealing or a genetic algorithm (see below)

Genetic Algorithms
Introduction

These are a powerful optimization technique, which rely on natural selection. There are many variants, and this is an area of active research. We will be holding a mini- workshop on these methods with the CEDC in February (keep an eye on our seminar page for details comp_model/seminars.html).

Algorithm

Set up a population of individuals, by encoding the variables in the problem as a set of genes.

Evaluate the fitness of the individuals, according to the objective function.

Breed the solutions, using a genetic operator.

Perform some random mutations in the gene pool.

Repeat from (2).

Again constraints may be incorporated by a repair algorithm after breeding / mutation (to ensure that the solutions remain feasible), or by using a penalty function.

Advantages of genetic algorithms

They require no knowledge or gradient information about the objective function, and they can handle discontinuities.

They can give information about the stability of a solution.

Can be used to give solutions for complex multi-dimensional optimization, for which few other techniques can give any information.

Disadvantages of genetic algorithms

They do not necessarily find the exact global optimum.

Require large number of evaluations of the objective function

A vanilla GA may not perform very well, it is necessary to incorporate problem specific information for optimal performance.

Summary

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

Certain problems are intrinsically "hard", but well understood: you should recognise them and choose the appropriate solution technique.

The techniques discussed and the problems which arise are exemplars for other areas of numerical analysis.

Last updated 3-Dec-97. Maintained by SJ Cox