Forthcoming events in this series
Clustering, reordering and random graphs
Abstract
From the point of view of a numerical analyst, I will describe some algorithms for:
- clustering data points based on pairwise similarity,
- reordering a sparse matrix to reduce envelope, two-sum or bandwidth,
- reordering nodes in a range-dependent random graph to reflect the range-dependency,
and point out some connections between seemingly disparate solution techniques. These datamining problems arise across a range of disciplines. I will mention a particularly new and important application from bioinformatics concerning the analysis of gene or protein interaction data.
Immersed interface methods for fluid dynamics problems
Abstract
Immersed interface methods have been developed for a variety of
differential equations on domains containing interfaces or irregular
boundaries. The goal is to use a uniform Cartesian grid (or other fixed
grid on simple domain) and to allow other boundaries or interfaces to
cut through this grid. Special finite difference formulas are developed
at grid points near an interface that incorporate the appropriate jump
conditions across the interface so that uniform second-order accuracy
(or higher) can be obtained. For fluid flow problems with an immersed
deformable elastic membrane, the jump conditions result from a balance
between the singular force imposed by the membrane, inertial forces if
the membrane has mass, and the jump in pressure across the membrane.
A second-order accurate method of this type for Stokes flow was developed
with Zhilin Li and more recently extended to the full incompressible
Navier-Stokes equations in work with Long Lee.
Inverse eigenvalue problems for quadratic matrix polynomials
Abstract
Feedback design for a second order control system leads to an
eigenstructure assignment problem for a quadratic matrix polynomial. It is
desirable that the feedback controller not only assigns specified
eigenvalues to the second order closed loop system, but also that the
system is robust, or insensitive to perturbations. We derive here new
sensitivity measures, or condition numbers, for the eigenvalues of the
quadratic matrix polynomial and define a measure of robustness of the
corresponding system. We then show that the robustness of the quadratic
inverse eigenvalue problem can be achieved by solving a generalized linear
eigenvalue assignment problem subject to structured perturbations.
Numerically reliable methods for solving the structured generalized linear
problem are developed that take advantage of the special properties of the
system in order to minimize the computational work required.
Modelling bilevel games in electricity
Abstract
Electricity markets facilitate pricing and delivery of wholesale power.
Generators submit bids to an Independent System Operator (ISO) to indicate
how much power they can produce depending on price. The ISO takes these bids
with demand forecasts and minimizes the total cost of power production
subject to feasibility of distribution in the electrical network.
\\
\\
Each generator can optimise its bid using a bilevel program or
mathematical program with equilibrium (or complementarity) constraints, by
taking the ISOs problem, which contains all generators bid information, at
the lower level. This leads immediately to a game between generators, where
a Nash equilibrium - at which each generator's bid maximises its profit
provided that none of the other generators changes its bid - is sought.
\\
\\
In particular, we examine the idealised model of Berry et al (Utility
Policy 8, 1999), which gives a bilevel game that can be modelled as an
"equilibrium problem with complementarity constraints" or EPCC.
Unfortunately, like bilevel games, EPCCs on networks may not have Nash
equilibria in the (common) case when one or more of links of the network is
saturated (at maximum capacity). Nevertheless we explore some theory and
algorithms for this problem, and discuss the economic implications of
numerical examples where equilibria are found for small electricity
networks.
Combinatorial structures in nonlinear programming
Abstract
Traditional optimisation theory and -methods on the basis of the
Lagrangian function do not apply to objective or constraint functions
which are defined by means of a combinatorial selection structure. Such
selection structures can be explicit, for example in the case of "min",
"max" or "if" statements in function evaluations, or implicit as in the
case of inverse optimisation problems where the combinatorial structure is
induced by the possible selections of active constraints. The resulting
optimisation problems are typically neither convex nor smooth and do not
fit into the standard framework of nonlinear optimisation. Users typically
treat these problems either through a mixed-integer reformulation, which
drastically reduces the size of tractable problems, or by employing
nonsmooth optimisation methods, such as bundle methods, which are
typically based on convex models and therefore only allow for weak
convergence results. In this talk we argue that the classical Lagrangian
theory and SQP methodology can be extended to a fairly general class of
nonlinear programs with combinatorial constraints. The paper is available
Exact real arithmetic
Abstract
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.
Generalised finite and infinite elements for flow acoustics
Improving spectral methods with optimized rational interpolation
Abstract
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.
Numerical issues arising in dynamic optimisation of process modelling applications
Abstract
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.
Eigenmodes of polygonal drums
Abstract
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.
Convergence analysis of linear and adjoint approximations with shocks
Geometry, PDEs fluid dynamics, and image processing
Abstract
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.)
Computing solutions of Laplace's equation by conformal mapping
Special Alan Curtis event
Abstract
- 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
On the convergence of interior point methods for linear programming
Abstract
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.
Spectral effects with quaternions
Abstract
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.
Computation of period orbits for the Navier-Stokes equations
Abstract
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.
On the solution of moving boundary value problems adaptive moving meshes
Superlinear convergence of conjugate gradients
Abstract
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.
Sobolev index estimation for hp-adaptive finite element methods
Abstract
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.
Recent results on accuracy and stability of numerical algorithms
Abstract
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.
Real symmetric matrices with multiple eigenvalues
Abstract
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.
Some complexity considerations in sparse LU factorization
Abstract
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.