Past Computational Mathematics and Applications Seminar

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.