Past Computational Mathematics and Applications Seminar

We calculate an analytic value for the correlation coefficient between a geometric, or exponential, Brownian motion and its time-average, a novelty being our use of divided differences to elucidate formulae. This provides a simple approximation for the value of certain Asian options regarding them as exchange options. We also illustrate that the higher moments of the time-average can be expressed neatly as divided differences of the exponential function via the Hermite-Genocchi integral relation, as well as demonstrating that these expressions agree with those obtained by Oshanin and Yor when the drift term vanishes.

The formation of steady patterns in one space dimension is generically

governed, at small amplitude, by the Ginzburg-Landau equation.

But in systems with a conserved quantity, there is a large-scale neutral

mode that must be included in the asymptotic analysis for pattern

formation near onset. The usual Ginzburg-Landau equation for the amplitude

of the pattern is then coupled to an equation for the large-scale mode.

\\

These amplitude equations show that for certain parameters all regular

periodic patterns are unstable. Beyond the stability boundary, there

exist stable stationary solutions in the form of spatially modulated

patterns or localised patterns. Many more exotic localised states are

found for patterns in two dimensions.

\\

Applications of the theory include convection in a magnetic field,

providing an understanding of localised states seen in numerical

simulations.

The lattice Boltzmann equation has been used successfully used to simulate

nearly incompressible flows using an isothermal equation of state, but

much less work has been done to determine stable implementations for other

equations of state. The commonly used nine velocity lattice Boltzmann

equation supports three non-hydrodynamic or "ghost'' modes in addition to

the macroscopic density, momentum, and stress modes. The equilibrium value

of one non-hydrodynamic mode is not constrained by the continuum equations

at Navier-Stokes order in the Chapman-Enskog expansion. Instead, we show

that it must be chosen to eliminate a high wavenumber instability. For

general barotropic equations of state the resulting stable equilibria do

not coincide with a truncated expansion in Hermite polynomials, and need

not be positive or even sign-definite as one would expect from arguments

based on entropy extremisation. An alternative approach tries to suppress

the instability by enhancing the damping the non-hydrodynamic modes using

a collision operator with multiple relaxation times instead of the common

single relaxation time BGK collision operator. However, the resulting

scheme fails to converge to the correct incompressible limit if the

non-hydrodynamic relaxation times are fixed in lattice units. Instead we

show that they must scale with the Mach number in the same way as the

stress relaxation time.

Consider approximating a set of discretely defined values $f_{1}, \ldots , f_{m}$ say at $x=x_{1}, x_{2}, \ldots, x_{m}$, with a chosen approximating form. Given prior knowledge that noise is present and that some might be outliers, a standard least squares approach based on $l_{2}$ norm of the error $\epsilon$ may well provide poor estimates. We instead consider a least squares approach based on a modified measure of the form $\tilde{\epsilon} = \epsilon (1+c^{2}\epsilon^{2})^{-\frac{1}{2}}$, where $c$ is a constant to be fixed.

\\

The choice of the constant $c$ in this estimator has a significant effect on the performance of the estimator both in terms of its algorithmic convergence to a solution and its ability to cope effectively with outliers. Given a prior estimate of the likely standard deviation of the noise in the data, we wish to determine a value of $c$ such that the estimator behaves like a robust estimator when outliers are present but like a least squares estimator otherwise.

\\

We describe approaches to determining suitable values of $c$ and illustrate their effectiveness on approximation with polynomial and radial basis functions. We also describe algorithms for computing the estimates based on an iteratively weighted linear least squares scheme.

We consider matrix groups defined in terms of scalar products. Examples of interest include the groups of

• complex orthogonal,
• real, complex, and conjugate symplectic,
• real perplectic,
• real and complex pseudo-orthogonal,
• pseudo-unitary

matrices. We

• Construct a variety of transformations belonging to these groups that imitate the actions of Givens rotations, Householder reflectors, and Gauss transformations.
• Describe applications for these structured transformations, including to generating random matrices in the groups.
• Show how to exploit group structure when computing the polar decomposition, the matrix sign function and the matrix square root on these matrix groups.

This talk is based on recent joint work with N. Mackey, D. S. Mackey, and N. J. Higham.

The heart can be described as an electrically driven mechanical pump. This

pump couldn't adapt to beat-by-beat changes in circulatory demand if there

was no feedback from the mechanical environment to the electrical control

processes. Cardiac mechano-electric feedback has been studied at various

levels of functional integration, from stretch-activated ion channels,

through mechanically induced changes in cardiac cells and tissue, to

clinically relevant observations in man, where mechanical stimulation of the

heart may either disturb or reinstate cardiac rhythmicity. The seminar will

illustrate the patho-physiological relevance of cardiac mechano-electric

feedback, introduce underlying mechanisms, and show the utility of iterating

between experimental research and mathematical modelling in studying this

phenomenon.

Semidefinite programming (SDP) techniques have been extremely successful

in many practical engineering design questions. In several of these

applications, the problem structure is invariant under the action of

some symmetry group, and this property is naturally inherited by the

underlying optimization. A natural question, therefore, is how to

exploit this information for faster, better conditioned, and more

reliable algorithms. To this effect, we study the associative algebra

associated with a given SDP, and show the striking advantages of a

careful use of symmetries. The results are motivated and illustrated

through applications of SDP and sum of squares techniques from networked

control theory, analysis and design of Markov chains, and quantum

information theory.

Discontinuous Galerkin methods (DG) use trial and test functions that are continuous within

elements and discontinuous at element boundaries. Although DG methods have been invented

in the early 1970s they have become very popular only recently.

\\

DG methods are very attractive for flow and transport problems in porous media since they

can be used to solve hyperbolic as well as elliptic/parabolic problems, (potentially) offer

high-order convergence combined with local mass balance and can be applied to unstructured,

non-matching grids.

\\

In this talk we present a discontinuous Galerkin method based on the non-symmetric interior

penalty formulation introduced by Wheeler and Rivi\`{e}re for an elliptic equation coupled to

a nonlinear parabolic/hyperbolic equation. The equations cover models for groundwater flow and

solute transport as well as two-phase flow in porous media.

\\

We show that the method is comparable in efficiency with the mixed finite element method for

elliptic problems with discontinuous coefficients. In the case of two-phase flow the method

can outperform standard finite volume schemes by a factor of ten for a five-spot problem and

also for problems with dominating capillary pressure.

The standard airframe industry tool for flutter analysis is based

on linear potential predictions of the aerodynamics. Despite the

limitations of the modelling this is even true in the transonic

range. There has been a heavy research effort in the past decade to

use CFD to generate the aerodynamics for flutter simulations, to

improve the reliability of predictions and thereby reduce the risk

and cost of flight testing. The first part of the talk will describe

efforts at Glasgow to couple CFD with structural codes to produce

a time domain simulation and an example calculation will be described for

the BAE SYSTEMS Hawk aircraft.

\\

\\

A drawback with time domain simulations is that unsteady CFD is still

costly and parametric searches to determine stability through the

growth or decay of responses can quickly become impractical. This has

motivated another active research effort in developing ways of

encapsulating the CFD level aerodynamic predictions in models which

are more affordable for routine application. A number of these

approaches are being developed (eg POD, system identification...)

but none have as yet reached maturity. At Glasgow effort has been

put into developing a method based on the behaviour of the

eigenspectrum of the discrete operator Jacobian, using Hopf

Bifurcation conditions to formulate an augmented system of

steady state equations which can be used to calculate flutter speeds

directly. The talk will give the first three dimensional example

of such a calculation.

\\

\\

For background reports on these topics see

http://www.aero.gla.ac.uk/Research/CFD/projects/aeroelastics/pubs/menu…

It is known for elliptic problems with smooth coefficients

that the solution is smooth in the interior of the domain;

low regularity is only possible near the boundary.

The $hp$-version of the FEM allows us to exploit this

property if we use meshes where the element size grows

porportionally to the element's distance to the boundary

and the approximation order is suitably linked to the

element size. In this way most degrees of freedom are

concentrated near the boundary.

\\

In this talk, we will discuss convergence and complexity

issues of the boundary concentrated FEM. We will show

that it is comparable to the classical boundary element

method (BEM) in that it leads to the same convergence rate

(error versus degrees of freedom). Additionally, it

generalizes the classical FEM since it does not require

explicit knowledge of the fundamental solution so that

it is also applicable to problems with (smooth) variable

coefficients.

Following the framework of Formaggia and Perotto (Numer.

Math. 2001 and 2003), anisotropic a posteriori error estimates have been

proposed for various elliptic and parabolic problems. The error in the

energy norm is bounded above by an error indicator involving the matrix

of the error gradient, the constant being independent of the mesh aspect

ratio. The matrix of the error gradient is approached using

Zienkiewicz-Zhu error estimator. Numerical experiments show that the

error indicator is sharp. An adaptive finite element algorithm which

aims at producing successive triangulations with high aspect ratio is

proposed. Numerical results will be presented on various problems such

as diffusion-convection, Stokes problem, dendritic growth.

While the average settling velocity of particles in a suspension has been successfully predicted, we are still unsuccessful with the r.m.s velocity, with theories suggesting a divergence with the size of

the container and experiments finding no such dependence. A possible resolution involves stratification originating from the spreading of the front between the clear liquid above and the suspension below. One theory describes the spreading front by a nonlinear diffusion equation

$\frac{\partial \phi}{\partial t} = D \frac{\partial }{\partial z}(\phi^{4/5}(\frac{\partial \phi}{\partial z})^{2/5})$.

\\

\\

Experiments and computer simulations find differently.

We consider the problem of recovering the position of a scattering surface

from measurements of the scattered field on a finite line above the surface.

A point source algorithm is proposed, based on earlier work by Potthast,

which reconstructs, in the first instance, the scattered field in the whole

region above the scattering surface. This information is used in a second

stage to locate the scatterer. We summarise the theoretical results that can

be obtained (error bounds on the reconstructed field as a function of the

noise level in the original measurements). For the case of a point source of

the incident field we present numerical experiments for both a steady source

(time harmonic excitation) and a pulse source typical of an antenna in

ground penetrating radar applications.

\\

This is joint work with Claire Lines (Brunel University).

The seminar will address issues related to numerical simulation

of non-linear behaviour of solid materials to impact loading.

The kinematic and constitutive aspects of the transition from

continuum to discontinuum will be presented as utilised

within an explicit finite element development framework.

Material softening, mesh sensitivity and regularisation of

solutions will be discussed.

To numerical analysts and other applied mathematicians Jacobians and Hessians

are matrices, i.e. rectangular arrays of numbers or algebraic expressions.

Possibly taking account of their sparsity such arrays are frequently passed

into library routines for performing various computational tasks.

\\

\\

A central goal of an activity called automatic differentiation has been the

accumulation of all nonzero entries from elementary partial derivatives

according to some variant of the chainrule. The elementary partials arise

in the user-supplied procedure for evaluating the underlying vector- or

scalar-valued function at a given argument.

\\

\\

We observe here that in this process a certain kind of structure that we

call "Jacobian scarcity" might be lost. This loss will make the subsequent

calculation of Jacobian vector-products unnecessarily expensive.

Instead we advocate the representation of the Jacobian as a linear computational

graph of minimal complexity. Many theoretical and practical questions remain unresolved.

Condition numbers are a central concept in numerical analysis.

They provide a natural parameter for studying the behavior of

algorithms, as well as sensitivity and geometric properties of a problem.

The condition number of a problem instance is usually a measure

of the distance to the set of ill-posed instances. For instance, the

classical Eckart and Young identity characterizes the condition

number of a square matrix as the reciprocal of its relative distance

to singularity.

\\

\\

We present concepts of conditioning for optimization problems and

for more general variational problems. We show that the Eckart and

Young identity has natural extension to much wider contexts. We also

discuss conditioning under the presence of block-structure, such as

that determined by a sparsity pattern. The latter has interesting

connections with the mu-number in robust control and with the sign-real

spectral radius.

The seminar will focus on mathematical modelling of physics phenomena,

with applications in e.g. mass and heat transfer, fluid flow, and

electromagnetic wave propagation. Simultaneous solutions of several

physics phenomena described by PDEs - multiphysics - will also be

presented and discussed.

\\

\\

All models will be realised through the use of the MATLAB based finite

element package FEMLAB. FEMLAB is a multiphysics modelling environment

with built-in PDE solvers for linear, non-linear, time dependent and

eigenvalue problems. For ease-of-use, it comprises ready-to-use applications

for various physics phenomena, and tailored applications for Structural

Mechanics, Electromagnetics, and Chemical Engineering. But in addition,

FEMLAB facilitates straightforward implementation of arbitrary coupled

non-linear PDEs, which brings about a great deal of flexibility in problem

definition. Please see http://ww.uk.comsol.com for more info.

\\

\\

FEMLAB is developed by the COMSOL Group, a Swedish headquartered spin-off

from the Royal Institute of Technology (KTH) in Stockholm, with offices

around the world. Its UK office is situated in The Oxford Science Park.

The process industries are one of the UK's major sectors and include

petrochemicals, pharmaceuticals, water, energy and the food industry,

amongst others. The design of a processing plant is a difficult task. This

is due to both the need to cater for multiple criteria (such as economics,

environmental and safety) and the use highly complex nonlinear models to

describe the behaviour of individual unit operations in the process. Early

in the design stages, an engineer may wish to use automated design tools to

generate conceptual plant designs which have potentially positive attributes

with respect to the main criteria. Such automated tools typically rely on

optimization for solving large mixed integer nonlinear programming models.

\\

\\

This talk presents an overview of some of the work done in the Computer

Aided Process Engineering group at UCL. Primary emphasis will be given to

recent developments in hybrid optimization methods, including the use of

graphical interfaces based on problem specific visualization techniques to

allow the engineer to interact with embedded optimization procedures. Case

studies from petrochemical and water industries will be presented to

demonstrate the complexities involved and illustrate the potential benefits

of hybrid approaches.

The simulation of sedimentary basins aims at reconstructing its historical

evolution in order to provide quantitative predictions about phenomena

leading to hydrocarbon accumulations. The kernel of this simulation is the

numerical solution of a complex system of time dependent, three

dimensional partial differential equations of mixed parabolic-hyperbolic

type in highly heterogeneous media. A discretisation and linearisation of

this system leads to large ill-conditioned non-symmetric linear systems

with three unknowns per mesh element.

\\

\\

In the seminar I will look at different preconditioning approaches for

these systems and at their parallelisation. The most effective

preconditioner which we developed so far consists in three stages: (i) a

local decoupling of the equations which (in addition) aims at

concentrating the elliptic part of the system in the "pressure block'';

(ii) an efficient preconditioning of the pressure block using AMG; (iii)

the "recoupling'' of the equations. Numerical results on real case

studies, exhibit (i) a significant reduction of sequential CPU times, up

to a factor 5 with respect to the current ILU(0) preconditioner, (ii)

robustness with respect to physical and numerical parameters, and (iii) a

speedup of up to 4 on 8 processors.

Rapidly oscillating problems, whether differential equations or

integrals, ubiquitous in applications, are allegedly difficult to

compute. In this talk we will endeavour to persuade the audience that

this is false: high oscillation, properly understood, is good for

computation! We describe methods for differential equations, based on

Magnus and Neumann expansions of modified systems, whose efficacy

improves in the presence of high oscillation. Likewise, we analyse

generalised Filon quadrature methods, showing that also their error

sharply decreases as the oscillation becomes more rapid.

In many applications of interest, such as the conformational

dynamics of molecules, large deterministic systems can exhibit

stochastic behaviour in a relative small number of coarse-grained

variables. This kind of dimension reduction, from a large deterministic

system to a smaller stochastic one, can be very useful in understanding

the problem. Whilst the subject of statistical mechanics provides

a wealth of explicit examples where stochastic models for coarse

variables can be found analytically, it is frequently the case

that applications of interest are not amenable to analytic

dimension reduction. It is hence of interest to pursue computational

algorithms for such dimension reduction. This talk will be devoted

to describing recent work on parameter estimation aimed at

problems arising in this context.

\\

\\

Joint work with Raz Kupferman (Jerusalem) and Petter Wiberg (Warwick)

A new filter method will be presented that attempts to find a feasible

point for sets of nonlinear sets of equalities and inequalities. The

method is intended to work for problems where the number of variables

or the number of (in)equalities is large, or both. No assumption is

made about convexity. The technique used is that of maintaining a list

of multidimensional "filter entries", a recent development of ideas

introduced by Fletcher and Leyffer. The method will be described, as

well as large scale numerical experiments with the corresponding

Fortran 90 module, FILTRANE.

In recent times, research into scattering of electromagnetic waves by complex objects

has assumed great importance due to its relevance to radar applications, where the

main objective is to identify targeted objects. In designing stealth weapon systems

such as military aircraft, control of their radar cross section is of paramount

importance. Aircraft in combat situations are threatened by enemy missiles. One

countermeasure which is used to reduce this threat is to minimise the radar cross

section. On the other hand, there is a demand for the enhancement of the radar cross

section of civilian spacecraft. Operators of communication satellites often request

a complicated differential radar cross section in order to assist with the tracking

of the satellite. To control the radar cross section, an essential requirement is a

capability for accurate prediction of electromagnetic scattering from complex objects.

\\

\\

One difficulty which is encountered in the development of suitable numerical solution

schemes is the existence of constraints which are in excess of those needed for a unique

solution. Rather than attempt to include the constraint in the equation set, the novel

approach which is presented here involves the use of the finite element method and the

construction of a specialised element in which the relevant solution variables are

appropriately constrained by the nature of their interpolation functions. For many

years, such an idea was claimed to be impossible. While the idea is not without its

difficulties, its advantages far outweigh its disadvantages. The presenter has

successfully developed such an element for primitive variable solutions to viscous

incompressible flows and wishes to extend the concept to electromagnetic scattering

problems.

\\

\\

Dr Mack has first degrees in mathematics and aeronautical engineering, plus a Masters

and a Doctorate, both in computational fluid dynamics. He has some thirty years

experience in this latter field. He pioneered the development of the innovative

solenoidal approach for the finite element solution of viscous incompressible flows.

At the time, such a radical idea was claimed in the literature to be impossible.

Much of this early research was undertaken during a six month sabbatical with the

Numerical Analysis Group at the Oxford University Computing Laboratory. Dr Mack has

since received funding from British Aerospace and the United States Department of

Defense to continue this research.

In addition to the announced topic of Pascal Matrices (abstract below) we will speak briefly about more recent work by Per-Olof Persson on generating simplicial meshes on regions defined by a function that gives the distance from the boundary. Our first goal was a short MATLAB code and we just submitted "A Simple Mesh Generator in MATLAB" to SIAM.

This is joint work with Alan Edelman at MIT and a little bit with Pascal. They had all the ideas.

Put the famous Pascal triangle into a matrix. It could go into a lower triangular L or its transpose L' or a symmetric matrix S:

These binomial numbers come from a recursion, or from the formula for i choose j, or functionally from taking powers of (1 + x).

The amazing thing is that L times L' equals S. (OK for 4 by 4) It follows that S has determinant 1. The matrices have other unexpected properties too, that give beautiful examples in teaching linear algebra. The proof of L L' = S comes 3 ways, I don't know which you will prefer:

1. By induction using the recursion formula for the matrix entries.
2. By an identity for the coefficients i+j choose j in S.
3. By applying both sides to the column vector [ 1 x x² x³ ... ]'.

The third way also gives a proof that S³ = -I but we doubt that result.

The rows of the "hypercube matrix" L² count corners and edges and faces and ... in n dimensional cubes.

Seminar moved to Week 8, 19 June 2003.