Past events

Is it possible to construct a computational model of the real numbers in which the sign

of every computed result is corrected determined? The answer is yes, both in theory and in

practice. The resulting viewpoint contrasts strongly with the traditional floating

point model. I will review the theoretical background and software design issues,

discuss previous attempts at implementation and finally demonstrate my own python and

C++ codes.

The pseudospectral method for solving boundary value problems on the interval

consists in replacing the solution by an interpolating polynomial in Lagrangian

form between well-chosen points and collocating at those same points.

\\

\\

Due to its globality, the method cannot handle steep gradients well (Markov's inequality).

We will present and discuss two means of improving upon this: the attachment of poles to

the ansatz polynomial, on one hand, and conformal point shifts on the other hand, both

optimally adapted to the problem to be solved.

Dynamic optimisation is a tool that enables the process industries to

compute optimal control strategies for important chemical processes.

Aspen DynamicsTM is a well-established commercial engineering software

package containing a dynamic optimisation tool. Its intuitive graphical

user interface and library of robust dynamic models enables engineers to

quickly and easily define a dynamic optimisation problem including

objectives, control vector parameterisations and constraints. However,

this is only one part of the story. The combination of dynamics and

non-linear optimisation can create a problem that can be very difficult

to solve due to a number of reasons, including non-linearities, poor

initial guesses, discontinuities and accuracy and speed of dynamic

integration. In this talk I will begin with an introduction to process

modelling and outline the algorithms and techniques used in dynamic

optimisation. I will move on to discuss the numerical issues that can

give us so much trouble in practice and outline some solutions we have

created to overcome some of them.

Many questions of interest to both mathematicians and physicists relate

to the behavior of eigenvalues and eigenmodes of the Laplace operator

on a polygon. Algorithmic improvements have revived the old "method

of fundamental solutions" associated with Fox, Henrici and Moler; is it

going to end up competitive with the state-of-the-art method of Descloux,

Tolley and Driscoll? This talk will outline the numerical issues but

give equal attention to applications including "can you hear the shape

of a drum?", localization of eigenmodes, eigenvalue avoidance, and the

design of drums that play chords.

\\

This is very much work in progress -- with graduate student Timo Betcke.

Image processing is an area with many important applications, as well as challenging problems for mathematicians. In particular, Fourier/wavelets analysis and stochastic/statistical methods have had major impact in this area. Recently, there has been increased interest in a new and complementary approach, using partial differential equations (PDEs) and differential-geometric models. It offers a more systematic treatment of geometric features of mages, such as shapes, contours and curvatures, etc., as well as allowing the wealth of techniques developed for PDEs and Computational Fluid Dynamics (CFD) to be brought to bear on image processing tasks.

I'll use two examples from my recent work to illustrate this synergy:

1. A unified image restoration model using Total Variation (TV) which can be used to model denoising, deblurring, as well as image inpainting (e.g. restoring old scratched photos). The TV idea can be traced to shock capturing methods in CFD and was first used in image processing by Rudin, Osher and Fatemi.

2. An "active contour" model which uses a variational level set method for object detection in scalar and vector-valued images. It can detect objects not necessarily defined by sharp edges, as well as objects undetectable in each channel of a vector-valued image or in the combined intensity. The contour can go through topological changes, and the model is robust to noise. The level set method was originally developed by Osher and Sethian for tracking interfaces in CFD.

(The above are joint works with Jackie Shen at the Univ. of Minnesota and Luminita Vese in the Math Dept at UCLA.)

2.00 pm Professor Iain Duff (RAL) Opening remarks

2.15 pm Professor M J D Powell (University of Cambridge)

Some developments of work with Alan on cubic splines

3.00 pm Professor Kevin Burrage (University of Queensland)

Stochastic models and simulations for chemically reacting systems

3.30 pm Tea/Coffee

4.00 pm Professor John Reid (RAL)

Sparse matrix research at Harwell and the Rutherford Appleton Laboratory

4.30 pm Dr Ian Jones (AEA PLC)

Computational fluid dynamics and the role of stiff solvers

5.00 pm Dr Lawrence Daniels (Hyprotech UK Ltd)

Current work with Alan on ODE solvers for HSL

Long-step primal-dual path-following algorithms constitute the

framework of practical interior point methods for

solving linear programming problems. We consider

such an algorithm and a second order variant of it.

We address the problem of the convergence of

the sequences of iterates generated by the two algorithms

to the analytic centre of the optimal primal-dual set.

Several real Lie and Jordan algebras, along with their associated

automorphism groups, can be elegantly expressed in the quaternion tensor

algebra. The resulting insight into structured matrices leads to a class

of simple Jacobi algorithms for the corresponding $n \times n$ structured

eigenproblems. These algorithms have many desirable properties, including

parallelizability, ease of implementation, and strong stability.

A method for computing periodic orbits for the Navier-Stokes

equations will be presented. The method uses a finite-element Galerkin

discretisation for the spatial part of the problem and a spectral

Galerkin method for the temporal part of the problem. The method will

be illustrated by calculations of the periodic flow behind a circular

cylinder in a channel. The problem has a simple reflectional symmetry

and it will be explained how this can be exploited to reduce the cost

of the computations.

The convergence of Krylov subspace methods like conjugate gradients

depends on the eigenvalues of the underlying matrix. In many cases

the exact location of the eigenvalues is unknown, but one has some

information about the distribution of eigenvalues in an asymptotic

sense. This could be the case for linear systems arising from a

discretization of a PDE. The asymptotic behavior then takes place

when the meshsize tends to zero.

\\

\\

We discuss two possible approaches to study the convergence of

conjugate gradients based on such information.

The first approach is based on a straightforward idea to estimate

the condition number. This method is illustrated by means of a

comparison of preconditioning techniques.

The second approach takes into account the full asymptotic

spectrum. It gives a bound on the asymptotic convergence factor

which explains the superlinear convergence observed in many situations.

This method is mathematically more involved since it deals with

potential theory. I will explain the basic ideas.

We develop an algorithm for estimating the local Sobolev regularity index

of a given function by monitoring the decay rate of its Legendre expansion

coefficients. On the basis of these local regularities, we design and

implement an hp--adaptive finite element method based on employing

discontinuous piecewise polynomials, for the approximation of nonlinear

systems of hyperbolic conservation laws. The performance of the proposed

adaptive strategy is demonstrated numerically.

The study of the finite precision behaviour of numerical algorithms dates back at least as far as Turing and Wilkinson in the 1940s. At the start of the 21st century, this area of research is still very active.

We focus on some topics of current interest, describing recent developments and trends and pointing out future research directions. The talk will be accessible to those who are not specialists in numerical analysis.

Specific topics intended to be addressed include

• Floating point arithmetic: correctly rounded elementary functions, and the fused multiply-add operation.
• The use of extra precision for key parts of a computation: iterative refinement in fixed and mixed precision.
• Gaussian elimination with rook pivoting and new error bounds for Gaussian elimination.
• Automatic error analysis.
• Application and analysis of hyperbolic transformations.

We describe "avoidance of crossing" and its explanation by von

Neumann and Wigner. We show Lax's criterion for degeneracy and then

discover matrices whose determinants give the discriminant of the

given matrix. This yields a simple proof of the bound given by

Ilyushechkin on the number of terms in the expansion of the discriminant

as a sum of squares. We discuss the 3 x 3 case in detail.

The talk will discuss unsymmetric sparse LU factorization based on

the Markowitz pivot selection criterium. The key question for the

author is the following: Is it possible to implement a sparse

factorization where the overhead is limited to a constant times

the actual numerical work? In other words, can the work be bounded

by o(sum(k, M(k)), where M(k) is the Markowitz count in pivot k.

The answer is probably NO, but how close can we get? We will give

several bad examples for traditional methods and suggest alternative

methods / data structure both for pivot selection and for the sparse

update operations.

We discuss two filters that are frequently used to smooth data.

One is the (nonlinear) median filter, that chooses the median

of the sample values in the sliding window. This deals effectively

with "outliers" that are beyond the correct sample range, and will

never be chosen as the median. A straightforward implementation of

the filter is expensive for large windows, particularly in two dimensions

(for images).

\\

\\

The second filter is linear, and known as "Savitzky-Golay". It is

frequently used in spectroscopy, to locate positions and peaks and

widths of spectral lines. This filter is based on a least-squares fit

of the samples in the sliding window to a polynomial of relatively

low degree. The filter coefficients are unlike the equiripple filter

that is optimal in the maximum norm, and the "maxflat" filters that

are central in wavelet constructions. Should they be better known....?

\\

\\

We will discuss the analysis and the implementation of both filters.

The asymptotic rate of convergence of the trapezium rule is

defined, for smooth functions, by the Euler-Maclaurin expansion.

The extension to other methods, such as Gauss rules, is straightforward;

this talk begins with some special cases, such as Periodic functions, and

functions with various singularities.

\\

\\

Convergence rates for ODEs (Initial and Boundary value problems)

and for PDEs are available, but not so well known. Extension to singular

problems seems to require methods specific to each situation. Some of

the results are unexpected - to me, anyway.

In the past few years we have developed some expertise in solving optimization

problems that involve large scale simulations in various areas of Computational

Geophysics and Engineering. We will discuss some of those applications here,

namely: inversion of seismic data, characterization of piezoelectrical crystals

material properties, optimal design of piezoelectrical transducers and

opto-electronic devices, and the optimal design of steel structures.

\\

\\

A common theme among these different applications is that the goal functional

is very expensive to evaluate, often, no derivatives are readily available, and

some times the dimensionality can be large.

\\

\\

Thus parallelism is a need, and when no derivatives are present, search type

methods have to be used for the optimization part. Additional difficulties can

be ill-conditioning and non-convexity, that leads to issues of global

optimization. Another area that has not been extensively explored in numerical

optimization and that is important in real applications is that of

multiobjective optimization.

\\

\\

As a result of these varied experiences we are currently designing a toolbox

to facilitate the rapid deployment of these techniques to other areas of

application with a minimum of retooling.

A-Posteriori Error estimates for high order Godunov finite

volume methods are presented which exploit the two solution

representations inherent in the method, viz. as piecewise

constants $u_0$ and cell-wise $q$-th order reconstructed

functions $R^0_q u_0$. The analysis provided here applies

directly to schemes such as MUSCL, TVD, UNO, ENO, WENO or any

other scheme that is a faithful extension of Godunov's method

to high order accuracy in a sense that will be made precise.

Using standard duality arguments, we construct exact error

representation formulas for derived functionals that are

tailored to the class of high order Godunov finite volume

methods with data reconstruction, $R^0_q u_0$. We then consider

computable error estimates that exploit the structure of higher

order Godunov finite volume methods. The analysis technique used

in this work exploits a certain relationship between higher

order Godunov methods and the discontinuous Galerkin method.

Issues such as the treatment of nonlinearity and the optional

post-processing of numerical dual data are also discussed.

Numerical results for linear and nonlinear scalar conservation

laws are presented to verify the analysis. Complete details can

be found in a paper appearing in the proceedings of FVCA3,

Porquerolles, France, June 24-28, 2002.

SMP (Symmetric Multi-Processors) hardware technologies are very popular

with vendors and end-users alike for a number of reasons. However, true

shared memory parallelism has experienced somewhat slower to take up

amongst the scientific-programming community. NAG has been at the

forefront of SMP technology for a number of years, and the NAG SMP

Library has shown the potential of SMP systems.

\\

\\

At the very high end, SMP hardware technologies are used as building

blocks of modern supercomputers, which truly consist of clusters of SMP

systems, for which no dedicated model of parallelism yet exists.

\\

\\

The aim of this talk is to introduce SMP systems and their potential.

Results from our work at NAG will also be introduced to show how SMP

parallelism, based on a shared memory paradigm, can be used to very

good effect and can produce high performance, scalable software. The

talk also aims to discuss some aspects of the apparent slow take up of

shared memory parallelism and the potential competition from PC (i.e.

Intel)-based cluster technology. The talk then aims to explore the

potential of SMP technology within "hybrid parallelism", i.e. mixed

distributed and shared memory modes, illustrating the point with some

preliminary work carried out by the author and others. Finally, a

number of potential future challenges to numerical analysts will be

discussed.

\\

\\

The talk is aimed at all who are interested in SMP technologies for

numerical computing, irrespective of any previous experience in the

field. The talk aims to stimulate discussion, by presenting some ideas,

backing these with data, not to stifle it in an ocean of detail!