Past Computational Mathematics and Applications Seminar

We present an O(N logN) algorithm for the calculation of the first N coefficients in an expansion of an analytic function in Legendre polynomials. In essence, the algorithm consists of an integration of a suitably weighted function along an ellipse, a task which can be accomplished with Fast Fourier Transform, followed by some post-processing.

This talk will focus on highly localized basis functions which exist for certain kernels and spaces associated with these kernels. Such kernels include certain radial basis functions (RBFs), their restrictions to spheres (SBFs), and their restrictions to more general manifolds embeddable in R^d. The first part of the talk will be of an introductory nature. It will discuss radial basis functions and their restriction to manifolds which give rise to various kernels on these manifolds. The talk will then focus on the development (for certain kernels) of highly localized Lagrange functions which serve as effective bases: i.e., bases which are stable and local. Scaled versions of these bases will then be used to establish the stability of the L² minimization operator in L^p, 1 ≤ p ≤ ∞, thus obtaining a multivariate analogue of a result of de Boor. Since these bases are scalable with the data, they have potential uses beyond approximation including meshless methods and, more generally, computations of a multiresolution nature. The talk is primarily based on joint work with T. Hangelbroek, F. J. Narcowich and X. Sun.

We discuss a class of high-order spectral-Galerkin surface integral algorithms with specific focus on simulating the scattering of electromagnetic waves by a collection of three dimensional deterministic and stochastic particles.

In the eighteenth century Gaspard Monge considered the problem of finding the best way of moving a pile of material from one site to another. This optimal transport problem has many applications such as mesh generation, moving mesh methods, image registration, image morphing, optical design, cartograms, probability theory, etc. The solution to an optimal transport problem can be found by solving the Monge-Amp\`{e}re equation, a highly nonlinear second order elliptic partial differential equation. Leonid Kantorovich, however, showed that it is possible to analyse optimal transport problems in a framework that naturally leads to a linear programming formulation. In recent years several efficient methods have been proposed for solving the Monge-Amp\`{e}re equation. For the linear programming problem, standard methods do not exploit the special properties of the solution and require a number of operations that is quadratic or even cubic in the number of points in the discretisation. In this talk I will discuss techniques that can be used to obtain more efficient methods.

Joint work with Chris Budd (University of Bath).

Meshless (or meshfree) methods are a relatively new numerical approach for the solution of ordinary- and partial differential equations. They offer the geometrical flexibility of finite elements but without requiring connectivity from the discretization support (ie a mesh). Meshless methods based on the collocation of radial basis functions (RBF methods) are particularly easy to code, and have a number of theoretical advantages as well as practical drawbacks. In this talk, an adaptive RBF scheme is presented for a novel application, namely the solution of (a rather broad class of) delayed- and neutral differential equations.

We consider saddle point problems arising as (linearized) optimality conditions in elliptic optimal control problems. The efficient solution of such systems is a core ingredient in second-order optimization algorithms. In the spirit of Bramble and Pasciak, the preconditioned systems are symmetric and positive definite with respect to a suitable scalar product. We extend previous work by Schoeberl and Zulehner and consider problems with control and state constraints. It stands out as a particular feature of this approach that an appropriate symmetric indefinite preconditioner can be constructed from standard preconditioners for those matrices which represent the inner products, such as multigrid cycles.

Numerical examples in 2D and 3D are given which illustrate the performance of the method, and limitations and open questions are addressed.

Interior-point methods for linear and convex quadratic programming

require the solution of a sequence of symmetric indefinite linear

systems which are used to derive search directions. Safeguards are

typically required in order to handle free variables or rank-deficient

Jacobians. We propose a consistent framework and accompanying

theoretical justification for regularizing these linear systems. Our

approach is akin to the proximal method of multipliers and can be

interpreted as a simultaneous proximal-point regularization of the

primal and dual problems. The regularization is termed "exact" to

emphasize that, although the problems are regularized, the algorithm

recovers a solution of the original problem. Numerical results will be

presented. If time permits we will illustrate current research on a

matrix-free implementation.

This is joint work with Michael Friedlander, University of British Columbia, Canada

Shifted Laplace preconditioners have attracted considerable attention as

a technique to speed up convergence of iterative solution methods for the

Helmholtz equation. In the talk we present a comprehensive spectral

analysis of the discrete 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 derive an optimal value for the shift, and to

explain the mesh-depency of the convergence of GMRES preconditioned

with a shifted Laplacian. We will illustrate our results with a seismic test

problem.

Joint work with: Yogi Erlangga (University of British Columbia) and Kees Vuik (TU Delft)

Nonlinear eigenvalue problem (NEP) is a class of eigenvalue problems where the matrix depends on the eigenvalue. We will first introduce some NEPs in real applications and some algorithms for general NEPs. Then we introduce our recent advances in NEPs, including second order Arnoldi algorithms for large scale quadratic eigenvalue problem (QEP), analysis and algorithms for symmetric eigenvalue problem with nonlinear rank-one updating, a new linearization for rational eigenvalue problem (REP).

This talk will introduce L1-regularized optimization problems that arise in image processing, and numerical methods for their solution. In particular, we will focus on methods of the split-Bregman type, which very efficiently solve large scale problems without regularization or time stepping. Applications include image

denoising, segmentation, non-local filters, and compressed sensing.

There is a widespread use of mathematical tools in finance and its

importance has grown over the last two decades. In this talk we

concentrate on optimization problems in finance, in particular on

numerical aspects. In this talk, we put emphasis on the mathematical problems and aspects, whereas all the applications are connected to the pricing of derivatives and are the

outcome of a cooperation with an international finance institution.

As one example, we take an in-depth look at the problem of hedging

barrier options. We review approaches from the literature and illustrate

advantages and shortcomings. Then we rephrase the problem as an

optimization problem and point out that it leads to a semi-infinite

programming problem. We give numerical results and put them in relation

to known results from other approaches. As an extension, we consider the

robustness of this approach, since it is known that the optimality is

lost, if the market data change too much. To avoid this effect, one can

formulate a robust version of the hedging problem, again by the use of

semi-infinite programming. The numerical results presented illustrate

the robustness of this approach and its advantages.

As a further aspect, we address the calibration of models being used in

finance through optimization. This may lead to PDE-constrained

optimization problems and their solution through SQP-type or

interior-point methods. An important issue in this context are

preconditioning techniques, like preconditioning of KKT systems, a very

active research area. Another aspect is the preconditioning aspect

through the use of implicit volatilities. We also take a look at the

numerical effects of non-smooth data for certain models in derivative

pricing. Finally, we discuss how to speed up the optimization for

calibration problems by using reduced order models.

Saddle-point problems occur frequently in liquid crystal modelling. For example, they arise whenever Lagrange multipliers are used for the pointwise-unit-vector constraints in director modelling, or in both general director and order tensor models when an electric field is present that stems from a constant voltage. Furthermore, in a director model with associated constraints and Lagrange multipliers, together with a coupled electric-field interaction, a particular ''double'' saddle-point structure arises. This talk will focus on a simple example of this type and discuss appropriate numerical solution schemes.

This is joint work with Eugene C. Gartland, Jr., Department of Mathematical Sciences, Kent State University.

We show that data assimilation using four-dimensional variation

(4DVar) can be interpreted as a form of Tikhonov regularisation, a

familiar method for solving ill-posed inverse problems. It is known from

image restoration problems that $L_1$-norm penalty regularisation recovers

sharp edges in the image better than the $L_2$-norm penalty

regularisation. We apply this idea to 4DVar for problems where shocks are

present and give some examples where the $L_1$-norm penalty approach

performs much better than the standard $L_2$-norm regularisation in 4DVar.

In dynamical systems given by an ODE, one is interested in the basin

of attraction of invariant sets, such as equilibria or periodic

orbits. The basin of attraction consists of solutions which converge

towards the invariant set. To determine the basin of attraction, one

can use a solution of a certain linear PDE which can be approximated

by meshless collocation.

The basin of attraction of an equilibrium can be determined through

sublevel sets of a Lyapunov function, i.e. a scalar-valued function

which is decreasing along solutions of the dynamical system. One

method to construct such a Lyapunov function is to solve a certain

linear PDE approximately using Meshless Collocation. Error estimates

ensure that the approximation is a Lyapunov function.

The basin of attraction of a periodic orbit can be analysed by Borg’s

criterion measuring the time evolution of the distance between

adjacent trajectories with respect to a certain Riemannian metric.

The sufficiency and necessity of this criterion will be discussed,

and methods how to compute a suitable Riemannian metric using

Meshless Collocation will be presented in this talk.

In the last few years, there has been renewed interest in stochastic

finite element methods (SFEMs), which facilitate the approximation

of statistics of solutions to PDEs with random data. SFEMs based on

sampling, such as stochastic collocation schemes, lead to decoupled

problems requiring only deterministic solvers. SFEMs based on

Galerkin approximation satisfy an optimality condition but require

the solution of a single linear system of equations that couples

deterministic and stochastic degrees of freedom. This is regarded as

a serious bottleneck in computations and the difficulty is even more

pronounced when we attempt to solve systems of PDEs with

random data via stochastic mixed FEMs.

In this talk, we give an overview of solution strategies for the

saddle-point systems that arise when the mixed form of the Darcy

flow problem, with correlated random coefficients, is discretised

via stochastic Galerkin and stochastic collocation techniques. For

the stochastic Galerkin approach, the systems are orders of

magnitude larger than those arising for deterministic problems. We

report on fast solvers and preconditioners based on multigrid, which

have proved successful for deterministic problems. In particular, we

examine their robustness with respect to the random diffusion

coefficients, which can be either a linear or non-linear function of

a finite set of random parameters with a prescribed probability

distribution.

In this talk we describe a general framework for deriving

modified equations for stochastic differential equations with respect to

weak convergence. We will start by quickly recapping of how to derive

modified equations in the case of ODEs and describe how these ideas can

be generalized in the case of SDEs. Results will be presented for first

order methods such as the Euler-Maruyama and the Milstein method. In the

case of linear SDEs, using the Gaussianity of the underlying solutions,

we will derive a SDE that the numerical method solves exactly in the

weak sense. Applications of modified equations in the numerical study

of Langevin equations and in the calculation of effective diffusivities

will also be discussed.

A brief introduction to the field of complex networks is carried out by means of some examples. Then, we focus on the topics of defining and applying centrality measures to characterise the nodes of complex networks. We combine this approach with methods for detecting communities as well as to identify good expansion properties on graphs. All these concepts are formally defined in the presentation. We introduce the subgraph centrality from a combinatorial point of view and then connect it with the theory of graph spectra. Continuing with this line we introduce some modifications to this measure by considering some known matrix functions, e.g., psi matrix functions, as well as new ones introduced here. Finally, we illustrate some examples of applications in particular the identification of essential proteins in proteomic maps.

Future discovery and control in biology and medicine will come from

the mathematical modeling of large-scale molecular biological data,

such as DNA microarray data, just as Kepler discovered the laws of

planetary motion by using mathematics to describe trends in

astronomical data. In this talk, I will demonstrate that

mathematical modeling of DNA microarray data can be used to correctly

predict previously unknown mechanisms that govern the activities of

DNA and RNA.

First, I will describe the computational prediction of a mechanism of

regulation, by using the pseudoinverse projection and a higher-order

singular value decomposition to uncover a genome-wide pattern of

correlation between DNA replication initiation and RNA expression

during the cell cycle. Then, I will describe the recent

experimental verification of this computational prediction, by

analyzing global expression in synchronized cultures of yeast under

conditions that prevent DNA replication initiation without delaying

cell cycle progression. Finally, I will describe the use of the

singular value decomposition to uncover "asymmetric Hermite functions,"

a generalization of the eigenfunctions of the quantum harmonic

oscillator, in genome-wide mRNA lengths distribution data.

These patterns might be explained by a previously undiscovered asymmetry

in RNA gel electrophoresis band broadening and hint at two competing

evolutionary forces that determine the lengths of gene transcripts.

Regularization techniques based on the Golub-Kahan iterative bidiagonalization belong among popular approaches for solving large discrete ill-posed problems. First, the original problem is projected onto a lower dimensional subspace using the bidiagonalization algorithm, which by itself represents a form of regularization by projection. The projected problem, however, inherits a part of the ill-posedness of the original problem, and therefore some form of inner regularization must be applied. Stopping criteria for the whole process are then based on the regularization of the projected (small) problem.

We consider an ill-posed problem with a noisy right-hand side (observation vector), where the noise level is unknown. We show how the information from the Golub-Kahan iterative bidiagonalization can be used for estimating the noise level. Such information can be useful for constructing efficient stopping criteria in solving ill-posed problems.

This is joint work by Iveta Hn\v{e}tynkov\'{a}, Martin Ple\v{s}inger, and Zden\v{e}k Strako\v{s} (Faculty of Mathematics and Physics, Charles University, and Institute of Computer Science, Academy of Sciences, Prague)

We consider rational approximations to the Faddeeva or plasma dispersion function, defined

as

$w(z) = e^{-z^{2}} \mbox{erfc} (-iz)$.

With many important applications in physics, good software for

computing the function reliably everywhere in the complex plane is required. In this talk

we shall derive rational approximations to $w(z)$ via quadrature, M\"{o}bius transformations, and best approximation. The various approximations are compared with regard to speed of convergence, numerical stability, and ease of generation of the coefficients of the formula.

In addition, we give preference to methods for which a single expression yields uniformly

high accuracy in the entire complex plane, as well as being able to reproduce exactly the

asymptotic behaviour

$w(z) \sim i/(\sqrt{\pi} z), z \rightarrow \infty$

(in an appropriate sector).

This is Joint work with: Stephan Gessner, St\'efan van der Walt

Invariant subspaces are a well-established tool in the theory of linear eigenvalue problems. They are also computationally more stable objects than single eigenvectors if one is interested in a group of closely clustered eigenvalues. A generalization of invariant subspaces to matrix polynomials can be given by using invariant pairs.

We investigate some basic properties of invariant pairs and give perturbation results, which show that invariant pairs have similarly favorable properties for matrix polynomials than do invariant subspaces have for linear eigenvalue problems. In the second part of the talk we discuss computational aspects, namely how to extract invariant pairs from linearizations of matrix polynomials and how to do efficient iterative refinement on them. Numerical examples are shown using the NLEVP collection of nonlinear eigenvalue test problems.

This talk is joint work with Daniel Kressner from ETH Zuerich.

Molecular Dynamics Simulations are a tool to study the behaviour

of atomic-scale systems. The simulations themselves solve the

equations of motion for hundreds to millions of particles over

thousands to billions of time steps. Due to the size of the

problems studied, such simulations are usually carried out on

large clusters or special-purpose hardware.

At a first glance, there is nothing much of interest for a

Numerical Analyst: the equations of motion are simple, the

integrators are of low order and the computational aspects seem

to focus on hardware or ever larger and faster computer

clusters.

The field, however, having been ploughed mainly by domain

scientists (e.g. Chemists, Biologists, Material Scientists) and

a few Computer Scientists, is a goldmine for interesting

computational problems which have been solved either badly or

not at all. These problems, although domain specific, require

sufficient mathematical and computational skill to make finding

a good solution potentially interesting for Numerical Analysts.

The proper solution of such problems can result in speed-ups

beyond what can be achieved by pushing the envelope on Moore's

Law.

In this talk I will present three examples where problems

interesting to Numerical Analysts arise. For the first two

problems, Constraint Resolution Algorithms and Interpolated

Potential Functions, I will present some of my own results. For

the third problem, using interpolations to efficiently compute

long-range potentials, I will only present some observations and

ideas, as this will be the main focus of my research in Oxford

and therefore no results are available yet.

CFD is an indispensible part of the design process for all major gas turbine components. The growth in the use of CFD from single-block structured mesh steady state solvers to highly resolved unstructured mesh unsteady solvers will be described, with examples of the design improvements that have been achieved. The European Commission has set stringent targets for the reduction of noise, emissions and fuel consumption to be achieved by 2020. The application of CFD to produce innovative designs to meet these targets will be described. The future direction of CFD towards whole engine simulations will also be discussed.