Past Computational Mathematics and Applications Seminar

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.

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.

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.

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.

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.

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.

This seminar has been cancelled.

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.

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.

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.

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.

No seminar today

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.

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.

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.

It has been known for some while now that every radial basis function

in common usage for multi-dimensional interpolation has associated with

it a uniquely defined Hilbert space, in which the radial basis function

is the `minimal norm interpolant'. This space is usually constructed by

utilising the positive definite nature of the radial function, but such

constructions turn out to be difficult to handle. We will present a

direct way of constructing the spaces, and show how to prove extension

theorems in such spaces. These extension theorems are at the heart of

improved error estimates in the $L_p$-norm.