Forthcoming events in this series
Asymptotic rates of convergence - for quadrature, ODEs and PDEs
Abstract
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.
A toolbox for optimal design
Abstract
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.
Analysis of some structured preconditioners for saddle point problems
A-Posteriori error estimates for higher order Godunov finite volume methods on unstructured meshes
Abstract
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 parallelism: Current achievements, future challenges
Abstract
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!
Computed tomography for X-rays: old 2-D results, relevance to new 3-D spiral CT problems
Oscillations in discrete solutions to the convection-diffusion equation
Abstract
It is well known that discrete solutions to the convection-diffusion
equation contain nonphysical oscillations when boundary layers are present
but not resolved by the discretisation. For the Galerkin finite element
method with linear elements on a uniform 1D grid, a precise statement as
to exactly when such oscillations occur can be made, namely, that for a
problem with mesh size h, constant advective velocity and different values
at the left and right boundaries, oscillations will occur if the mesh
P\'{e}clet number $P_e$ is greater than one. In 2D, however, the situation
is not so well understood. In this talk, we present an analysis of a 2D
model problem on a square domain with grid-aligned flow which enables us
to clarify precisely when oscillations occur, and what can be done to
prevent them. We prove the somewhat surprising result that there are
oscillations in the 2D problem even when $P_e$ is less than 1. Also, we show that there
are distinct effects arising from differences in the top and bottom
boundary conditions (equivalent to those seen in 1D), and the non-zero
boundaries parallel to the flow direction.
Algebraic modeling systems and mathematical programming
Abstract
Algebra based modeling systems are becoming essential elements in the
application of large and complex mathematical programs. These systems
enable the abstraction, expression and translation of practical
problems into reliable and effective operational systems. They provide
the bridged between algorithms and real world problems by automating
the problem analysis and translation into specific data structures and
provide computational services required by different solvers. The
modeling system GAMS will be used to illustrate the design goals and
main features of such systems. Applications in use and under
development will be used to provide the context for discussing the
changes in user focus and future requirements. This presents new sets
of opportunities and challenges to the supplier and implementer of
mathematical programming solvers and modeling systems.
Adaptive finite elements for optimal control
Abstract
A systematic approach to error control and mesh adaptation for
optimal control of systems governed by PDEs is presented.
Starting from a coarse mesh, the finite element spaces are successively
enriched in order to construct suitable discrete models.
This process is guided by an a posteriori error estimator which employs
sensitivity factors from the adjoint equation.
We consider different examples with the stationary Navier-Stokes
equations as state equation.
On the condition number of bases in Banach spaces
Iterative methods for PDE eigenvalue problems
Abstract
Finite Element approximation of surfactant spreading on a thin film
A new preconditioning technique for the solution of the biharmonic problem
Abstract
In this presentation we examine the convergence characteristics of a
Krylov subspace solver preconditioned by a new indefinite
constraint-type preconditioner, when applied to discrete systems
arising from low-order mixed finite element approximation of the
classical biharmonic problem. The preconditioning operator leads to
preconditioned systems having an eigenvalue distribution consisting of
a tightly clustered set together with a small number of outliers. We
compare the convergence characteristics of a new approach with the
convergence characteristics of a standard block-diagonal Schur
complement preconditioner that has proved to be extremely effective in
the context of mixed approximation methods.
\\
\\
In the second part of the presentation we are concerned with the
efficient parallel implementation of proposed algorithm on modern
shared memory architectures. We consider use of the efficient parallel
"black-box'' solvers for the Dirichlet Laplacian problems based on
sparse Cholesky factorisation and multigrid, and for this purpose we
use publicly available codes from the HSL library and MGNet collection.
We compare the performance of our algorithm with sparse direct solvers
from the HSL library and discuss some implementation related issues.
Distribution tails of condition numbers for the polyhedral conic feasibility problem
Abstract
(Joint work with Felipe Cucker and Dennis Cheung, City University of Hong Kong.)
\\
\\
Condition numbers are important complexity-theoretic tools to capture
a "distillation" of the input aspects of a computational problem that
determine the running time of algorithms for its solution and the
sensitivity of the computed output. The motivation for our work is the
desire to understand the average case behaviour of linear programming
algorithms for a large class of randomly generated input data in the
computational model of a machine that computes with real numbers. In
this model it is not known whether linear programming is polynomial
time solvable, or so-called "strongly polynomial". Closely related to
linear programming is the problem of either proving non-existence of
or finding an explicit example of a point in a polyhedral cone defined
in terms of certain input data. A natural condition number for this
computational problem was developed by Cheung and Cucker, and we analyse
its distributions under a rather general family of input distributions.
We distinguish random sampling of primal and dual constraints
respectively, two cases that necessitate completely different techniques
of analysis. We derive the exact exponents of the decay rates of the
distribution tails and prove various limit theorems of complexity
theoretic importance. An interesting result is that the existence of
the k-th moment of Cheung-Cucker's condition number depends only very
mildly on the distribution of the input data. Our results also form
the basis for a second paper in which we analyse the distributions of
Renegar's condition number for the randomly generated linear programming
problem.
Eigenvalues of Locally Perturbed Toeplitz Matrices
Abstract
Toeplitz matrices enjoy the dual virtues of ubiquity and beauty. We begin this talk by surveying some of the interesting spectral properties of such matrices, emphasizing the distinctions between infinite-dimensional Toeplitz matrices and the large-dimensional limit of the corresponding finite matrices. These basic results utilize the algebraic Toeplitz structure, but in many applications, one is forced to spoil this structure with some perturbations (e.g., by imposing boundary conditions upon a finite difference discretization of an initial-boundary value problem). How do such
perturbations affect the eigenvalues? This talk will address this question for "localized" perturbations, by which we mean perturbations that are restricted to a single entry, or a block of entries whose size remains fixed as the matrix dimension grows. One finds, for a broad class of matrices, that sufficiently small perturbations fail to alter the spectrum, though the spectrum is exponentially sensitive to other perturbations. For larger real single-entry
perturbations, one observes the perturbed eigenvalues trace out curves in the complex plane. We'll show a number of illustrations of this phenomenon for tridiagonal Toeplitz matrices.
\\
\\
This talk describes collaborative work with Albrecht Boettcher, Marko Lindner, and Viatcheslav Sokolov of TU Chemnitz.
Solution of massive support vector machine problems
Abstract
We investigate the use of interior-point and semismooth methods for solving
quadratic programming problems with a small number of linear constraints,
where the quadratic term consists of a low-rank update to a positive
semi-definite matrix. Several formulations of the support vector machine
fit into this category. An interesting feature of these particular problems
is the volume of data, which can lead to quadratic programs with between 10
and 100 million variables and, if written explicitly, a dense $Q$ matrix.
Our codes are based on OOQP, an object-oriented interior-point code, with the
linear algebra specialized for the support vector machine application.
For the targeted massive problems, all of the data is stored out of core and
we overlap computation and I/O to reduce overhead. Results are reported for
several linear support vector machine formulations demonstrating that the
methods are reliable and scalable and comparing the two approaches.
Spectral element methods for viscoelastic flow problems
Spectral inclusion and spectral exactness for non-selfadjoint differential equation eigenproblems
Abstract
Non-selfadjoint singular differential equation eigenproblems arise in a number of contexts, including scattering theory, the study of quantum-mechanical resonances, and hydrodynamic and magnetohydrodynamic stability theory.
\\
\\
It is well known that the spectra of non-selfadjoint operators can be pathologically sensitive to perturbation of the
operator. Wilkinson provides matrix examples in his classical text, while Trefethen has studied the phenomenon extensively through pseudospectra, which he argues are often of more physical relevance than the spectrum itself. E.B. Davies has studied the phenomenon particularly in the context of Sturm-Liouville operators and has shown that the eigenfunctions and associated functions of non-selfadjoint singular Sturm-Liouville operators may not even form a complete set in $L^2$.
\\
\\
In this work we ask the question: under what conditions can one expect the regularization process used for selfadjoint singular Sturm-Liouville operators to be successful for non-selfadjoint operators? The answer turns out to depend in part on the so-called Sims Classification of the problem. For Sims Case I the process is not guaranteed to work, and indeed Davies has very recently described the way in which spurious eigenvalues may be generated and converge to certain curves in the complex plane.
\\
\\
Using the Titchmarsh-Weyl theory we develop a very simple numerical procedure which can be used a-posteriori to distinguish genuine eigenvalues from spurious ones. Numerical results indicate that it is able to detect not only the spurious eigenvalues due to the regularization process, but also spurious eigenvalues due to the numerics on an already-regular problem. We present applications to quantum mechanical resonance calculations and to the Orr-Sommerfeld equation.
\\
\\
This work, in collaboration with B.M. Brown in Cardiff, has recently been generalized to Hamiltonian systems.
Numerical methods for stiff systems of ODEs
Abstract
Stiff systems of ODEs arise commonly when solving PDEs by spectral methods,
so conventional explicit time-stepping methods require very small time steps.
The stiffness arises predominantly through the linear terms, and these
terms can be handled implicitly or exactly, permitting larger time steps.
This work develops and investigates a class of methods known as
'exponential time differencing'. These methods are shown to have a
number of advantages over the more well-known linearly implicit
methods and integrating factor methods.
The Kestrel interface to the NEOS server
Abstract
The Kestrel interface for submitting optimization problems to the NEOS Server augments the established e-mail, socket, and web interfaces by enabling easy usage of remote solvers from a local modeling environment.
\\
\\
Problem generation, including the run-time detection of syntax errors, occurs on the local machine using any available modeling language facilities. Finding a solution to the problem takes place on a remote machine, with the result returned in the native modeling language format for further processing. A byproduct of the Kestrel interface is the ability to solve multiple problems generated by a modeling language in parallel.
\\
\\
This mechanism is used, for example, in the GAMS/AMPL solver available through the NEOS Server, which internally translates a submitted GAMS problem into AMPL. The resulting AMPL problem is then solved through the NEOS Server via the Kestrel interface. An advantage of this design is that the GAMS to AMPL translator does not need to be collocated with the AMPL solver used, removing restrictions on solver choice and reducing administrative costs.
\\
\\
This talk is joint work with Elizabeth Dolan.
Tridiagonal matrices and trees
Abstract
Tridiagonal matrices and three term recurrences and second order equations appear amazingly often, throughout all of mathematics. We won't try to review this subject. Instead we look in two less familiar directions.
\\
\\
Here is a tridiagonal matrix problem that waited surprisingly long for a solution. Forward elimination factors T into LDU, with the pivots in D as usual. Backward elimination, from row n to row 1, factors T into U_D_L_. Parlett asked for a proof that diag(D + D_) = diag(T) + diag(T^-1).^-1. In an excellent paper (Lin Alg Appl 1997) Dhillon and Parlett extended this four-diagonal identity to block tridiagonal matrices, and also applied it to their "Holy Grail" algorithm for the eigenproblem. I would like to make a different connection, to the Kalman filter.
\\
\\
The second topic is a generalization of tridiagonal to "tree-diagonal". Unlike the interval, the tree can branch. In the matrix T, each vertex is connected only to its neighbors (but a branch point has more than two neighbors). The continuous analogue is a second order differential equation on a tree. The "non-jump" conditions at a meeting of N edges are continuity of the potential (N-1 equations) and Kirchhoff's Current Law (1 equation). Several important properties of tridiagonal matrices, including O(N) algorithms, survive on trees.