Past Computational Mathematics and Applications Seminar

Radial basis functions have been used for decades for the interpolation of scattered,

high-dimensional data. Recently they have attracted interest as methods for simulating

partial differential equations as well. RBFs do not require a grid or triangulation, they

offer the possibility of spectral accuracy with local refinement, and their implementation

is very straightforward. A number of theoretical and practical breakthroughs in recent years

has improved our understanding and application of these methods, and they are currently being

tested on real-world applications in shallow water flow on the sphere and tear film evolution

in the human eye.

Joint work with Yogi Erlangga and Kees Vuik.

Shifted Laplace preconditioners have attracted considerable attention as a technique to speed up convergence of iterative solution methods for the Helmholtz equation. In this paper we present a comprehensive spectral analysis of the Helmholtz operator preconditioned with a shifted Laplacian. Our analysis is valid under general conditions. The propagating medium can be heterogeneous, and the analysis also holds for different types of damping, including a radiation condition for the boundary of the computational domain. By combining the results of the spectral analysis of the preconditioned Helmholtz operator with an upper bound on the GMRES-residual norm we are able to provide an optimal value for the shift, and to explain the mesh-depency of the convergence of GMRES preconditioned with a shifted Laplacian. We illustrate our results with a seismic test problem.

Many large-scale optimization problems arise in the context of the discretization of infinite dimensional applications. In such cases, the description of the finite-dimensional problem is not unique, but depends on the discretization used, resulting in a natural multi-level description. How can such a problem structure be exploited, in discretized problems or more generally? The talk will focus on discussing this issue in the context of unconstrained optimization and in relation with the classical multigrid approach to elliptic systems of partial differential equations. Both theoretical convergence properties of special purpose algorithms and their numerical performances will be discussed. Perspectives will also be given.

Collaboration with S. Gratton, A. Sartenaer and M. Weber.

The talk starts with a general introduction of the convex

quadratic semidefinite programming problem (QSDP), followed by a number of

interesting examples arising from finance, statistics and computer sciences.

We then discuss the concept of primal nondegeneracy for QSDP and show that

some QSDPs are nondegenerate and others are not even under the linear

independence assumption. The talk then moves on to the algorithmic side by

introducing the dual approach and how it naturally leads to Newton's method,

which is quadratically convergent for degenerate problems. On the

implementation side of the Newton method, we stress that direct methods for

the linear equations in Newton's method are impossible simply because the

equations are quite large in size and dense. Our numerical experiments use

the conjugate gradient method, which works quite well for the nearest

correlation matrix problem. We also discuss difficulties for us to find

appropriate preconditioners for the linear system encountered. The talk

concludes in discussing some other available methods and some future topics.

The aim of this talk is to describe several methods for numerically approximating

the integral of a multivariate highly oscillatory function. We begin with a review

of the asymptotic and Filon-type methods developed by Iserles and Nørsett. Using a

method developed by Levin as a point of departure we will construct a new method that

uses the same information as the Filon-type method, and obtains the same asymptotic

order, while not requiring moments. This allows us to integrate over nonsimplicial

domains, and with complicated oscillators.

High-performance computing is an important tool for computational science.

Oxford University has recently decided to invest £3M in a new supercomputing

facility which is under development now. In this seminar I will give an overview

of some activities in Oxford and provide a vision for supercomputing here.

I will discuss some of the numerical analysis software and tools,

such as Distributed Matlab and indicate some of the challenges at

the intersection of numerical analysis and high-performance computing.

The aim of this talk is to give some understanding of the theory of matrix $p$'th roots (solutions to the nonlinear matrix equation $X^{p} = A$), to explain how and how not to compute roots, and to describe some applications. In particular, an application in finance will be described concerning roots of transition matrices from Markov models.

One of the oldest approach in meshless methods for PDEs is the Interpolating Moving Least Squares (IMLS) technique developed in the 1980s. Although widely accepted by users working in fields as diverse as geoinformatics and crack dynamics I shall take a fresh look at this method and ask for the equivalent difference operators which are generated implicitly. As it turns out, these operators are optimal only in trivial cases and are "strange" in general. I shall try to exploit two different approaches for the computation of these operators.

On the other hand (and very different from IMLS), Total Variation Flow (TVF) PDEs are the most recent developments in image processing and have received much attention lately. Again I shall show that they are able to generate "strange" discrete operators and that they easily can behave badly although they may be properly implemented.

We will explore connections between the structure of high-dimensional convex polytopes and information acquisition for compressible signals. A classical result in the field of convex polytopes is that if N points are distributed Gaussian iid at random in dimension n<<N, then only order (log N)^n of the points are vertices of their convex hull. Recent results show that provided n grows slowly with N, then with high probability all of the points are vertices of its convex hull. More surprisingly, a rich "neighborliness" structure emerges in the faces of the convex hull. One implication of this phenomenon is that an N-vector with k non-zeros can be recovered computationally efficiently from only n random projections with n=2e k log(N/n). Alternatively, the best k-term approximation of a signal in any basis can be recovered from 2e k log(N/n) non-adaptive measurements, which is within a log factor of the optimal rate achievable for adaptive sampling. Additional implications for randomized error correcting codes will be presented.

This work was joint with David L. Donoho.

In this talk we present different strategies for regularization of the pure Newton method

(minimization problems)and of the Gauss-Newton method (systems of nonlinear equations).

For these schemes, we prove general convergence results. We establish also the global and

local worst-case complexity bounds. It is shown that the corresponding search directions can

be computed by a standard linear algebra technique.

We present a novel enhanced finite element method for the Darcy problem starting from the non stable

continuous $P_1 / P_0$ finite element spaces enriched with multiscale functions. The method is a departure

from the standard mixed method framework used in these applications. The methods are derived in a Petrov-Galerkin

framework where both velocity and pressure trial spaces are enriched with functions based on residuals of strong

equations in each element and edge partition. The strategy leads to enhanced velocity space with an element of

the lowest order Raviart-Thomas space and to a stable weak formulation preserving local mass conservation.

Numerical tests validate the method.

Jointly with Gabriel R Barrenechea, Universidad de Concepcion &amp;

Frederic G C Valentin, LNCC

Strong horizontal gradients of density are responsible for the occurence of a large number of (often catastrophic) flows, such as katabatic winds, dust storms, pyroclastic flows and powder-snow avalanches. For a large number of applications, the overall density contrast in the flow remains small and simulations are carried in the Boussinesq limit, where density variations only appear in the body-force term. However, pyroclastic flows and powder-snow avalanches involve much larger density contrasts, which implies that the inhomogeneous Navier-Stokes equations need to be solved, along with a closure equation describing the mass diffusion. We propose a Lagrange-Galerkin numerical scheme to solve this system, and prove optimal error bounds subject to constraints on the order of the discretization and the time-stepping. Simulations of physical relevance are then shown.

An integral part of the brain is a fluid flow system that is separate from brain tissue and the cerebral blood flow system: cerebrospinal fluid (CSF) is produced near the centre of the brain, flows out and around the brain, including around the spinal cord and is absorbed primarily in a region between the brain tissue and the skull. Hydrocephalus covers a broad range of anomalous flow and pressure situations: the normal flow path can become blocked, other problems can occur which result in abnormal tissue deformation or pressure changes. This talk will describe work that treats brain tissue as a poroelastic matrix through which the CSF can move when normal flow paths are blocked, producing tissue deformation and pressure changes. We have a number of models, the simplest treating the brain and CSF flow as having spherial symmetry ranging to more complex, fully three-dimensional computations. As well as considering acute hydrocephalus, we touch on normal pressure hydrocephalus, idiopathic intracranial hypertension and simulation of an infusion test. The numerical methods used are a combination of finite difference and finite element techniques applied to an interesting set of hydro-elastic equations.

In this talk, we report on recent activities in the development of automatic differentiation tools for Matlab and CapeML, a common intermediate language for process control, and highlight some recent AD applications. Lastly, we show the potential for parallelisation created by AD and comment on the impact on scientific computing due to emerging multicore chips which are providing substantial thread-based parallelism in a "pizza box" form factor.

We consider the solution of sequences of linear systems by preconditioned iterative methods. Such systems arise, for example, in applications such as CFD and structural mechanics. In some cases it is important to avoid the recomputation of preconditioners for subsequent systems. We propose an algebraic strategy that replaces new preconditioners by old preconditioners with simple updates. Efficiency of the new strategy, which generalizes the approach of Benzi and Bertaccini, is demonstrated using numerical experiments.

This talk presents results of joint work with Jurjen Duintjer Tebbens.

The use of discrete differential forms in the construction of finite element discretisations of the Sobolev spaces H^s, H(div) and H(curl) is now routinely applied by numerical analysts and engineers alike. However, little attention has been paid to the conditioning of the resulting stiffness matrices, particularly in the case of the non-uniform meshes that arise when adaptive refinement algorithms are used. We study this issue and show that the matrices are generally rather poorly conditioned. Typically, diagonal scaling is applied (often unwittingly) as a preconditioner. However, whereas diagonal scaling removes the effect of the mesh non-uniformity in the case of Sobolev spaces H^s, we show this is not so in the case of the spaces H(curl) and H(div). We trace the reason behind this difference, and give a simple remedy for curing the problem.

We propose a new parallel domain decomposition algorithm to solve symmetric linear systems of equations derived from the discretization of PDEs on general unstructured grids of triangles or tetrahedra. The algorithm is based on a single-level Schwarz alternating procedure and a modified conjugate gradient solver. A single layer of overlap has been adopted in order to simplify the data-structure and minimize the overhead. This approach makes the global convergence rate vary slightly with the number of domains and the algorithm becomes highly scalable. The algorithm has been implemented in Fortran 90 using MPI and hence portable to different architectures. Numerical experiments have been carried out on a SunFire 15K parallel computer and have been shown superlinear performance in some cases.

Non-normal matrices can be tiresome; some eigenvalues may be phlegmatic while others may be volatile. Computable error bounds are rarely used in such computations. We offer a way to proceed. Let (e,q,p') be an approximate eigentriple for non-normal B. Form column and row residuals r = Bq - qe and s' = p'B - ep'. We establish the relation between the smallest perturbation E, in both spectral and Frobenius norms, that makes the approximations correct and the norms of r and s'. Our results extend to the case when q and p are tall thin matrices and e is a small square matrix. Now regard B as a perturbation of B-E to obtain a (first order) bound on the error in e as a product of ||E|| and the condition number of e, namely (||q|| ||p'||)/|p'q|.

The use of preconditioned Newton-Krylov methods is in many applications mandatory for computing efficiently the solution of large nonlinear systems of equations. However, the available preconditioners are often sub-optimal, due to the changing nature of the linearized operator. This the case, for instance, for quasi-Newton methods where the Jacobian (and its preconditioner) are kept fixed at each non-linear iteration, with the rate of convergence usually degraded from quadratic to linear. Updated Jacobians, on the other hand require updated preconditioners, which may not be readily available. In this work we introduce an adaptive preconditioning technique based on the Krylov subspace information generated at previous steps in the nonlinear iteration. In particular, we use to advantage an adaptive technique suggested for restarted GMRES to enhance existing preconditioners with information about (almost) invariant subspaces constructed by GMRES at previous stages in the nonlinear iteration. We provide guidelines on the choice of invariant-subspace basis used in the construction of our preconditioner and demonstrate the improved performance on various test problems. As a useful general application we consider the case of augmented systems preconditioned by block triangular matrices based on the structure of the system matrix. We show that a sufficiently good solution involving the primal space operator allows for an efficient application of our adaptive technique restricted to the space of dual variables.

In this talk we will review classical multigrid methods and give an overview on algebraic multigrid methods, in particular the "classical" approach to AMG by Ruge and Stueben.

After that we will introduce a new class of multilevel methods. These new AMGs on one hand and exploit information based on filtering vectors and on the other hand, information about the inverse matrix is used to drive the coarsening process.

This new kind of AMG will be discussed and compared with "classical" AMG from a theoretical point of view as well as by showing some numerical examples.

We present a new stabilized mixed finite element method for the linear elasticity problem. The approach is based on the introduction of Galerkin least-squares terms arising from the constitutive and equilibrium equations, and from the relation defining the rotation in terms of the displacement.

We show that the resulting augmented variational formulation and the associated Galerkin scheme are well posed, and that the latter becomes locking-free and asymptotically locking-free for Dirichlet and mixed boundary conditions, respectively. In particular, the discrete scheme allows the utilization of Raviart-Thomas spaces of lowest order for the stress tensor, piecewise linear elements for the displacement, and piecewise constants for the rotation.

In the case of mixed boundary conditions, the essential one (Neumann) is imposed weakly, which yields the introduction of the trace of the displacement as a suitable Lagrange multiplier. This trace is then approximated by piecewise linear elements on an independent partition of the Neumann boundary whose mesh size needs to satisfy a compatibility condition with the mesh size associated to the triangulation of the domain. A reliable and efficient a-posteriori error estimate is also described. Finally, several numerical results illustrating the performance of the augmented scheme are reported.