Forthcoming events in this series
Spectral asymptotics of the damped wave operator: theory, simulations and applications
Augmented linear systems - methods and observations
Abstract
The talk will focus on solution methods for augmented linear systems of
the form
\\
\\
$[ A B ][x] = [b] [ B' 0 ][y] [0]$.
\\
\\
Augmented linear systems of this type arise in several areas of
numerical applied mathematics including mixed finite element / finite
difference discretisations of flow equations (Darcy flow and Stokes
flow), electrical network simulation and optimisation. The general
properties of such systems are that they are large, sparse and
symmetric, and efficient solution techniques should make use of the
block structure inherent in the system as well as of these properties.
\\
\\
Iterative linear solution methods will be described that
attempt to take advantage of the structure of the system, and
observations on augmented
systems, in particular the distribution of their eigenvalues, will be
presented which lead to further iterative methods and also to
preconditioners for existing solution methods. For the particular case
of Darcy flow, comments on properties of domain decomposition methods
of additive Schwarz type and similarities to incomplete factorisation
preconditioners will be made.
An adaptive finite element algorithm for the solution of time-dependent free-surface incompressible flow problems
Self-scaled barriers for semidefinite programming
Abstract
I am going to show that all self-scaled barriers for the
cone of symmetric positive semidefinite matrices are of the form
$X\mapsto -c_1\ln\det X +c_0$ for some constants $c_1$ > $0,c_0 \in$ \RN.
Equivalently one could state say that all such functions may be
obtained via a homothetic transformation of the universal barrier
functional for this cone. The result shows that there is a certain
degree of redundancy in the axiomatic theory of self-scaled barriers,
and hence that certain aspects of this theory can be simplified. All
relevant concepts will be defined. In particular I am going to give
a short introduction to the notion of self-concordance and the
intuitive ideas that motivate its definition.
An efficient Schur preconditioner based on modified discrete wavelet transforms
Exception-free arithmetic on the extended reals
Abstract
Interval arithmetic is a way to produce guaranteed enclosures of the
results of numerical calculations. Suppose $f(x)$ is a real
expression in real variables $x= (x_1, \ldots, x_n)$, built up from
the 4 basic arithmetic operations and other 'standard functions'. Let
$X_1, \ldots, X_n$ be (compact) real intervals. The process of {\em
interval evaluation} of $f(X_1, ..., X_n)$ replaces each real
operation by the corresponding interval operation wherever it occurs
in $f$, e.g. $A \times B$ is the smallest interval containing $\{a
\times b \mid a \in A, b \in B\}$.
\\
\\
As is well known, it yields a guaranteed enclosure for the true range
$\{f(x_1, \ldots, x_n) \mid x_1 \in X_1, \ldots, x_n \in X_n\}$,
provided no exceptions such as division by (an interval containing)
zero occur during evaluation.
\\
\\
Interval arithmetic takes set inputs and produces set outputs. Noting
this, we show there is a consistent way to extend arithmetic to $R^* =
R \cup \{-\infty, +\infty\}$, such that interval evaluation continues
to give enclosures, and there are {\em no exceptions}. The basic
ideas are: the usual set-theory meaning of evaluating a relation at a
set; and taking topological closure of the graph of a function in a
suitable $(R^{*})^n$. It gives rigorous meaning to intuitively
sensible statements like $1/0 = \{-\infty, +\infty\}$, $0/0 = R^*$
(but $(x/x)_{|x=0} = 1$), $\sin(+\infty) = [-1,1]$, etc.
\\
\\
A practical consequence is that an exception-free floating-point
interval arithmetic system is possible. Such a system is implemented
at hardware level in the new Sun Fortran compiler, currently on
beta-release.
Analysis of the Cholesky method with iterative refinement for solving the symmetric definite generalized eigenproblem
Abstract
The Cholesky factorization approach to solving the symmetric definite generalized eigenvalue problem
$Ax = \lambda Bx$, where $A$ is symmetric and $B$ is symmetric positive definite, computes a Cholesky factorization $B = LL^T$ and solves the equivalent standard symmetric eigenvalue problem $C y = \l y$ where $C = L^{-1} A L^{-T}$. Provided that a stable eigensolver is used, standard error analysis says that the computed eigenvalues are exact for $A+\dA$ and $B+\dB$ with $\max( \normt{\dA}/\normt{A}, \normt{\dB}/\normt{B} )$
bounded by a multiple of $\kappa_2(B)u$, where $u$ is the unit roundoff. We take the Jacobi method as the eigensolver and explain how backward error bounds potentially much smaller than $\kappa_2(B)u$ can be obtained.
To show the practical utility of our bounds we describe a vibration problem from structural engineering in which $B$ is ill conditioned yet the error bounds are small. We show how, in cases of instability, iterative refinement based on Newton's method can be used to produce eigenpairs with small backward errors.
Our analysis and experiments also give insight into the popular Cholesky--QR method used in LAPACK, in which the QR method is used as the eigensolver.
C*-algebras and pseudospectra of large Toeplitz matrices
Abstract
In contrast to spectra, pseudospectra of large Toeplitz matrices
behave as nicely as one could ever expect. We demonstrate some
basic phenomena of the asymptotic distribution of the spectra
and pseudospectra of Toeplitz matrices and show how by employing
a few simple $C^*$-algebra arguments one can prove rigorous
convergence results for the pseudospectra. The talk is a survey
of the development since a 1992 paper by Reichel and Trefethen
and is not addressed to specialists, but rather to a general
mathematically interested audience.
Sensitivity analysis for design and control in an elastic CAD-free framework for multi-model configurations
Abstract
This lecture is about the extension of our CAD-Free platform to the simulation and sensitivity analysis for design and control of multi-model configurations. We present the different ingredients of the platform using simple models for the physic of the problem. This presentation enables for an easy evaluation of coupling, design and control strategies.
\\
\\
Sensitivity analysis has been performed using automatic differentiation in direct or reverse modes and finite difference or complex variable methods. This former approach is interesting for cases where only one control parameter is involved as we can evaluate the state and sensitivity in real time with only one evaluation.
\\
\\
We show that for some class of applications, incomplete sensitivities can be evaluated dropping the state dependency which leads to a drastic cost reduction.
\\
\\
The design concerns the structural characteristic of the model and our control approach is based in perturbating the inflow incidence by a time dependent flap position described by a dynamic system encapsulating a gradient based minimization algorithm expressed as dynamic system. This approach enables for the definition of various control laws for different minimization algorithm.
\\
\\
We present dynamic minimization algorithms based on the coupling of several heavy ball method. This approach enables for global minimization at a cost proportional to the number of balls times the cost of one steepest descent minimization.
\\
\\
Two and three dimension flutter problem simulation and control are presented and the sensitivity of the different parameters in the model is discussed.
Continuation and bifurcation analysis of periodic solutions of partial differential equations
Abstract
There is a growing interest in the study of periodic phenomena in
large-scale nonlinear dynamical systems. Often the high-dimensional
system has only low-dimensional dynamics, e.g., many reaction-diffusion
systems or Navier-Stokes flows at low Reynolds number. We present an
approach that exploits this property in order to compute branches of
periodic solutions of the large system of ordinary differential
equations (ODEs) obtained after a space discretisation of the PDE. We
call our approach the Newton-Picard method. Our method is based on the
recursive projection method (RPM) of Shroff and Keller but extends this
method in many different ways. Our technique tries to combine the
performance of straightforward time integration with the advantages of
solving a nonlinear boundary value problem using Newton's method and a
direct solver. Time integration works well for very stable limit
cycles. Solving a boundary value problem is expensive, but works also
for unstable limit cycles.
\\
\\
We will present some background material on RPM. Next we will explain
the basic features of the Newton-Picard method for single shooting. The
linearised system is solved by a combination of direct and iterative
techniques. First, we isolate the low-dimensional subspace of unstable
and weakly stable modes (using orthogonal subspace iteration) and
project the linearised system on this subspace and on its
(high-dimensional) orthogonal complement. In the high-dimensional
subspace we use iterative techniques such as Picard iteration or GMRES.
In the low-dimensional (but "hard") subspace, direct methods such as
Gaussian elimination or a least-squares are used. While computing the
projectors, we also obtain good estimates for the dominant,
stability-determining Floquet multipliers. We will present a framework
that allows us to monitor and steer the convergence behaviour of the
method.
\\
\\
RPM and the Newton-Picard technique have been developed for PDEs that
reduce to large systems of ODEs after space discretisation. In fact,
both methods can be applied to any large system of ODEs. We will
indicate how these methods can be applied to the discretisation of the
Navier-Stokes equations for incompressible flow (which reduce to an
index-2 system of differential-algebraic equations after space
discretisation when written in terms of velocity and pressure.)
\\
\\
The Newton-Picard method has already been extended to the computation
of bifurcation points on paths of periodic solutions and to multiple
shooting. Extension to certain collocation and finite difference
techniques is also possible.
Improvements for iterative methods?
Abstract
Krylov subspace methods offer good possibilities for the solution of
large sparse linear systems of equations.For general systems, some of
the popular methods often show an irregular type of convergence
behavior and one may wonder whether that could be improved or not. Many
suggestions have been made for improvement and the question arises
whether these corrections are cosmetic or not. There is also the
question whether the irregularity shows inherent numerical instability.
In such cases one should take extra care in the application of
smoothing techniques. We will discuss strategies that work well and
strategies that might have been expected to work well.
Entropy Splitting for High-Order Numerical Simulation of Compressible Turbulence
Abstract
This work forms part of a larger research project to develop efficient
low-dissipative high-order numerical techniques for high-speed
turbulent flow simulation, including shock wave interactions with
turbulence. The requirements on a numerical method are stringent.For
the turbulence the method must be capable of resolving accurately a
wide range of length scales, whilst for shock waves the method must be
stable and not generate excessive local oscillations. Conventional
methods are either too dissipative, or incapable of shock capturing.
Higher-order ENO, WENO or hybrid schemes are too expensive for
practical computations. Previous work of Yee, Sandham & Djomehri
(1999) developed high-order shock-capturing schemes which minimize the
use of numerical dissipation away from shock
waves. The objective of the present study is to further minimize the
use of numerical dissipation for shock-free compressible turbulence
simulations.
Cheap Newton steps for discrete time optimal control problems: automatic differentiation and Pantoja's algorithm
Abstract
In 1983 Pantoja described a stagewise construction of the exact Newton
direction for a discrete time optimal control problem. His algorithm
requires the solution of linear equations with coefficients given by
recurrences involving second derivatives, for which accurate values are
therefore required.
\\
\\
Automatic differentiation is a set of techniques for obtaining derivatives
of functions which are calculated by a program, including loops and
subroutine calls, by transforming the text of the program.
\\
\\
In this talk we show how automatic differentiation can be used to
evaluate exactly the quantities required by Pantoja's algorithm,
thus avoiding the labour of forming and differentiating adjoint
equations by hand.
\\
\\
The cost of calculating the newton direction amounts to the cost of
solving one set of linear equations, of the order of the number of
control variables, for each time step. The working storage cost can be made
smaller than that required to hold the solution.
Analysis of a mean field model of superconducting vortices
Preconditioning constrained systems
Abstract
The general importance of preconditioning in combination with an
appropriate iterative technique for solving large scale linear(ised)
systems is widely appreciated. For definite problems (where the
eigenvalues lie in a half-plane) there are a number of preconditioning
techniques with a range of applicability, though there remain many
difficult problems. For indefinite systems (where there are eigenvalues
in both half-planes), techniques are generally not so well developed.
Constraints arise in many physical and mathematical problems and
invariably give rise to indefinite linear(ised) systems: the incompressible
Navier-Stokes equations describe conservation of momentum in the
presence of viscous dissipation subject to the constraint of
conservation of mass, for transmission problems the solution on an
interior domain is often solved subject to a boundary integral which
imposes the exterior field, in optimisation the appearance of
constraints is ubiquitous...
\\
\\
We will describe two approaches to preconditioning such constrained
systems and will present analysis and numerical results for each. In
particular, we will describe the applicability of these techniques to
approximations of incompressible Navier-Stokes problems using mixed
finite element approximation.
Arithmetic on the European Logarithmic Microprocessor
Abstract
As an alternative to floating-point, several workers have proposed the use
of a logarithmic number system, in which a real number is represented as a
fixed-point logarithm. Multiplication and division therefore proceed in
minimal time with no rounding error. However, the system can only offer an
overall advantage if addition and subtraction can be performed with speed
and accuracy at least equal to that of floating-point, but this has
hitherto been difficult to achieve. We will present a number of original
techniques by which this has now been accomplished. We will then
demonstrate by means of simulations that the logarithmic system offers
around twofold improvements in speed and accuracy, and finally will
describe a new European collaborative project which aims to develop a
logarithmic microprocessor during the next three years.
On the convergence of an implicitly restarted Arnoldi method
Abstract
We show that Sorensen's (1992) implicitly restarted Arnoldi method
(IRAM) (including its block extension) is non-stationary simultaneous
iteration in disguise. By using the geometric convergence theory for
non-stationary simultaneous iteration due to Watkins and Elsner (1991)
we prove that an implicitly restarted Arnoldi method can achieve a
super-linear rate of convergence to the dominant invariant subspace of
a matrix. We conclude with some numerical results the demonstrate the
efficiency of IRAM.