Past Computational Mathematics and Applications Seminar

Sparse grids yield numerical solutions to PDEs with a

significantly reduced number of degrees of freedom. The relative

benefit increases with the dimensionality of the problem, which makes

multi-factor models for financial derivatives computationally tractable.

An outline of a convergence proof for the so called combination

technique will be given for a finite difference discretisation of the

heat equation, for which sharp error bounds can be shown.

Numerical examples demonstrate that by an adaptive (heuristic)

choice of the subspaces European and American options with up to thirty

(and most likely many more) independent variables can be priced with

high accuracy.

The computation of flows of compressible fluids will be

discussed, exploiting the symmetric form of the equations describing

compressible flow.

A cell is a wonderously complex object. In this talk I will

give an overview of some of the mathematical frameworks that are needed

in order to make progress to understanding the complex dynamics of a

cell. The talk will consist of a directed random walk through discrete

Markov processes, stochastic differential equations, anomalous diffusion

and fractional differential equations.

An essential first step in many problems of numerical analysis and

computer graphics is to cover a region with a reasonably regular mesh.

We describe a short MATLAB code that begins with a "distance function"

to describe the region: $d(x)$ is the distance to the boundary

(with d

The talk will start with an introduction to recent development in Mathematica, with emphasis on numerical computing. This will be followed by a discussion of graph drawing algorithms for the display of relational information, in particular force directed algorithms. The talk will show that by employing multilevel approach and octree data structure, it is possible to achieve fast display of very large relational information, without compromising the quality.

We consider the solution of a linear system of equations using the GMRES iterative method. In some applications, performing inexact matrix-vector products in this method may be interesting, provided that a reasonable convergence of GMRES is achieved. A GMRES algorithm where the matrix vector product is performed inexactly is termed ”inexact GMRES algorithm”. An application of this idea occurs in computational electromagnetics, where the fast multipole method provides approximations of the matrix-vector product within a user-defined precision, and where these inaccurate matrix-vector products are all the more cheaper (in terms of CPU time) as the user-defined precision is low. The key point is then to design a control strategy of the accuracy of the matrix-vector product so that the GMRES converges better in the sense that 1) the inexact method achieves a satisfactory limiting accuracy, 2) within a reasonable number of steps.

\\

In [1], a relaxation strategy is proposed for general systems and validated on a large set of numerical experiments. This work is based on heuristic considerations and proposes a strategy that enables a convergence of the GMRES iterates $x_{k}$ within a relative normwise backward error $\frac{\|b−Ax_{k}\|}{\|A\| \|x_{k}\| + \|b\|}$ less than a prescribed quantity $\eta$ > 0, on a significant number of numerical experiments. Similar strategies have been applied to the solution of device simulation problems using domain decomposition [2] and to the preconditioning of a radiation diffusion problem in [5].

\\

A step toward a theoretical explanation of the observed behaviour of the inexact GMRES is proposed in [3, 4]. In this talk, we show that in spite of this considerable theoretical study, the experimental work of [1] is not fully understood yet. We give an overview of the questions that still remains open both in exact arithmetic and in floating-point arithmetic, and we provide some insights into the solution of some of them. Possible applications of this work for the preconditioned GMRES method, when the matrix-vector product is accurate but the preconditioning operation is approximate, are also investigated, based on [3].

\\

\\

References

\\

[1] A. Bouras and V. Frayss´e. Inexact matrix-vector products in Krylov methods for solving linear systems: a relaxation strategy. SIAM Journal on Matrix Analysis and Applications, 2004. To appear.

\\

[2] A. Bouras, V. Frayss´e, and L. Giraud. A relaxation strategy for inner-outer linear solvers in domain decomposition methods. Technical Report TR/PA/00/17, CERFACS, Toulouse, France, 2000.

\\

[3] V. Simoncini and D. B. Szyld. Theory of inexact Krylov subspace methods and applications to scientific computing. SIAM Journal Scientific Computing, 25:454–477, 2003.

\\

[4] J. van den Eshof and G. L. G. Sleijpen. Inexact Krylov subspace methods for linear systems. SIAM Journal on Matrix Analysis and Applications, February 2004. To appear.

\\

[5] J. S. Warsa, M. Benzi, T. A. Warein, and J. E. Morel. Preconditioning a mixed discontinuous finite element method for radiation diffusion. Numerical Linear Algebra with Applications, 2004. To appear.

In recent years, there has been considerable interest, especially in the context of

fluid-dynamics, in nonconforming finite element methods that are based on discontinuous

piecewise polynomial approximation spaces; such approaches are referred to as discontinuous

Galerkin (DG) methods. The main advantages of these methods lie in their conservation properties, their ability to treat a wide range of problems within the same unified framework, and their great flexibility in the mesh-design. Indeed, DG methods can easily handle non-matching grids and non-uniform, even anisotropic, polynomial approximation degrees. Moreover, orthogonal bases can easily be constructed which lead to diagonal mass matrices; this is particularly advantageous in unsteady problems. Finally, in combination with block-type preconditioners, DG methods can easily be parallelized.

\\

\\

In this talk, we introduce DG discretizations of mixed field and potential-based formulations of

eddy current problems in the time-harmonic regime. For the electric field formulation, the

divergence-free constraint within non-conductive regions is imposed by means of a Lagrange

multiplier. This allows for the correct capturing of edge and corner singularities in polyhedral domains; in contrast, additional Sobolev regularity must be assumed in the DG formulation, and their conforming counterparts, when regularization techniques are employed. In particular, we present a mixed method involving discontinuous $P^\ell-P^\ell$ elements, which includes a normal jump stabilization term, and a non-stabilized variant employing discontinuous $P^\ell-P^{\ell+1}$ elements.The first formulation delivers optimal convergence rates for the vector-valued unknowns in a suitable energy norm, while the second (non-stabilized) formulation is designed to yield optimal convergence rates in both the $L^2$--norm, as well as in a suitable energy norm. For this latter method, we also develop the {\em a posteriori} error estimation of the mixed DG approximation of the Maxwell operator. Indeed, by employing suitable Helmholtz decompositions of the error, together with the conservation properties of the underlying method, computable upper bounds on the error, measured in terms of the energy norm, are derived.

\\

\\

Numerical examples illustrating the performance of the proposed methods will be presented; here,

both conforming and non-conforming (irregular) meshes will be employed. Our theoretical and

numerical results indicate that the proposed DG methods provide promising alternatives to standard conforming schemes based on edge finite elements.

The satisfiability (SAT) problem is a central problem in mathematical

logic, computing theory, and artificial intelligence. We consider

instances of SAT specified by a set of boolean variables and a

propositional formula in conjunctive normal form. Given such an instance,

the SAT problem asks whether there is a truth assignment to the variables

such that the formula is satisfied. It is well known that SAT is in

general NP-complete, although several important special cases can be

solved in polynomial time. Extending the work of de Klerk, Warners and van

Maaren, we present new linearly sized semidefinite programming (SDP)

relaxations arising from a recently introduced paradigm of higher

semidefinite liftings for discrete optimization problems. These

relaxations yield truth assignments satisfying the SAT instance if a

feasible matrix of sufficiently low rank is computed. The sufficient rank

values differ between relaxations and can be viewed as a measure of the

relative strength of each relaxation. The SDP relaxations also have the

ability to prove that a given SAT formula is unsatisfiable. Computational

results on hard instances from the SAT Competition 2003 show that the SDP

approach has the potential to complement existing techniques for SAT.

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.