I will present an analysis of a continuous version of the compressed sensing problem, where the l^1 norm is replaced by the total variation of measures, and one aims to recover the positions and amplitudes of Dirac masses. We show that provided that the Diracs are sufficiently separated under a Fisher metric (which accounts for the geometry of the problem), stable recovery can be achieved when the number of random samples scales linearly with sparsity (up to log factors). This is joint work with Nicolas Keriven and Gabriel Peyre.

# Past Computational Mathematics and Applications Seminar

Iterative algorithms for the solution of convex quadratic optimization problems are investigated, which exploit inaccurate matrix-vector products. Theoretical bounds on the performance of a Conjugate Gradients method are derived, the necessary quantities occurring in the theoretical bounds estimated and a new practical algorithm derived. Numerical experiments suggest that the new method has significant potential, including in the steadily more important context of multi-precision computations.

Everybody is familiar with the concept of eigenvalues of a matrix. In this talk, we consider the nonlinear eigenvalue problem. These are problems for which the eigenvalue parameter appears in a nonlinear way in the equation. In physics, the Schroedinger equation for determining the bound states in a semiconductor device, introduces terms with square roots of different shifts of the eigenvalue. In mechanical and civil engineering, new materials often have nonlinear damping properties. For the vibration analysis of such materials, this leads to nonlinear functions of the eigenvalue in the system matrix.

One particular example is the sandwhich beam problem, where a layer of damping material is sandwhiched between two layers of steel. Another example is the stability analysis of the Helmholtz equation with a noise excitation produced by burners in a combustion chamber. The burners lead to a boundary condition with delay terms (exponentials of the eigenvalue).

We often receive the question: “How can we solve a nonlinear eigenvalue problem?” This talk explains the different steps to be taken for using Krylov methods. The general approach works as follows: 1) approximate the nonlinearity by a rational function; 2) rewrite this rational eigenvalue problem as a linear eigenvalue problem and then 3) solve this by a Krylov method. We explain each of the three steps.

The fully coupled numerical simulation of different physical processes, which can typically occur

at variable time and space scales, is often a very challenging task. A common feature of such models is that

their discretization gives rise to systems of linearized equations with an inherent block structure, which

reflects the properties of the set of governing PDEs. The efficient solution of a sequence of systems with

matrices in a block form is usually one of the most time- and memory-demanding issue in a coupled simulation.

This effort can be carried out by using either iteratively coupled schemes or monolithic approaches, which

tackle the problem of the system solution as a whole.

This talk aims to discuss recent advances in the monolithic solution of coupled multi-physics problems, with

application to poromechanical simulations in fractured porous media. The problem is addressed either by proper

sparse approximations of the Schur complements or special splittings that can partially uncouple the variables

related to different physical processes. The selected approaches can be included in a more general preconditioning

framework that can help accelerate the convergence of Krylov subspace solvers. The generalized preconditioner

relies on approximately decoupling the different processes, so as to address each single-physics problem

independently of the others. The objective is to provide an algebraic framework that can be employed as a

general ``black-box'' tool and can be regarded as a common starting point to be later specialized for the

particular multi-physics problem at hand.

Numerical experiments, taken from real-world examples of poromechanical problems and fractured media, are used to

investigate the behaviour and the performance of the proposed strategies.

Domain decomposition methods are widely employed for the numerical solution of partial differential equations on parallel computers. We develop an adjoint-based a posteriori error analysis for overlapping multiplicative Schwarz domain decomposition and for overlapping additive Schwarz. In both cases the numerical error in a user-specified functional of the solution (quantity of interest), is decomposed into a component that arises due to the spatial discretization and a component that results from of the finite iteration between the subdomains. The spatial discretization error can be further decomposed in to the errors arising on each subdomain. This decomposition of the total error can then be used as part of a two-stage approach to construct a solution strategy that efficiently reduces the error in the quantity of interest.

The development of effective solvers for high frequency wave propagation problems, such as those described by the Helmholtz equation, presents significant challenges. One promising class of solvers for such problems are parallel domain decomposition methods, however, an appropriate coarse space is typically required in order to obtain robust behaviour (scalable with respect to the number of domains, weakly dependant on the wave number but also on the heterogeneity of the physical parameters). In this talk we introduce a coarse space based on generalised eigenproblems in the overlap (GenEO) for the Helmholtz equation. Numerical results within FreeFEM demonstrate convergence that is effectively independent of the wave number and contrast in the heterogeneous coefficient as well as good performance for minimal overlap.

Can we get rigorous answers when computing with real and complex numbers? There are now many applications where this is possible thanks to a combination of tools from computer algebra and traditional numerical computing. I will give an overview of such methods in the context of two projects I'm developing. The first project, Arb, is a library for arbitrary-precision ball arithmetic, a form of interval arithmetic enabling numerical computations with rigorous error bounds. The second project, Fungrim, is a database of knowledge about mathematical functions represented in symbolic form. It is intended to function both as a traditional reference work and as a software library to support symbolic-numeric methods for problems involving transcendental functions. I will explain a few central algorithmic ideas and explain the research goals of these projects.

Partial differential equations (PDEs) are among the most universal tools used in modelling problems in nature and man-made complex systems. For example, stochastic PDEs are a fundamental ingredient in models for nonlinear filtering problems in chemical engineering and weather forecasting, deterministic Schroedinger PDEs describe the wave function in a quantum physical system, deterministic Hamiltonian-Jacobi-Bellman PDEs are employed in operations research to describe optimal control problems where companys aim to minimise their costs, and deterministic Black-Scholes-type PDEs are highly employed in portfolio optimization models as well as in state-of-the-art pricing and hedging models for financial derivatives. The PDEs appearing in such models are often high-dimensional as the number of dimensions, roughly speaking, corresponds to the number of all involved interacting substances, particles, resources, agents, or assets in the model. For instance, in the case of the above mentioned financial engineering models the dimensionality of the PDE often corresponds to the number of financial assets in the involved hedging portfolio. Such PDEs can typically not be solved explicitly and it is one of the most challenging tasks in applied mathematics to develop approximation algorithms which are able to approximatively compute solutions of high-dimensional PDEs. Nearly all approximation algorithms for PDEs in the literature suffer from the so-called "curse of dimensionality" in the sense that the number of required computational operations of the approximation algorithm to achieve a given approximation accuracy grows exponentially in the dimension of the considered PDE. With such algorithms it is impossible to approximatively compute solutions of high-dimensional PDEs even when the fastest currently available computers are used. In the case of linear parabolic PDEs and approximations at a fixed space-time point, the curse of dimensionality can be overcome by means of Monte Carlo approximation algorithms and the Feynman-Kac formula. In this talk we introduce new nonlinear Monte Carlo algorithms for high-dimensional nonlinear PDEs. We prove that such algorithms do indeed overcome the curse of dimensionality in the case of a general class of semilinear parabolic PDEs and we thereby prove, for the first time, that a general semilinear parabolic PDE with a nonlinearity depending on the PDE solution can be solved approximatively without the curse of dimensionality.

The Landau-DeGennes Q-model of uniaxial nematic liquid crystals seeks a rank-one

traceless tensor Q that minimizes a Frank-type energy plus a double well potential

that confines the eigenvalues of Q to lie between -1/2 and 1. We propose a finite

element method (FEM) which preserves this basic structure and satisfies a discrete

form of the fundamental energy estimates. We prove that the discrete problem Gamma

converges to the continuous one as the meshsize tends to zero, and propose a discrete

gradient flow to compute discrete minimizers. Numerical experiments confirm the ability

of the scheme to approximate configurations with half-integer defects, and to deal with

colloidal and electric field effects. This work, joint with J.P. Borthagaray and S.

Walker, builds on our previous work for the Ericksen's model which we review briefly.