material
Numerical Analysis and Scientific Computing Seminars 2008/09
Semester One

26 Sep
2008 Recent Progress in Computing the Matrix Exponential and its Frechet Derivative
Nick Higham (University of Manchester)
3.00  Frank Adams Room 1, Alan Turing BuildingAbstract (click to view)I will begin by introducing the concept of functions of matrices and briefly describing the history of the subject. Then I will describe the state of the art version of the most widely used method for computing $e^A$, the scaling and squaring method, and show how it can be extended to compute both $e^A$ and the Fr\'echet derivative at $A$ in the direction $E$ at a cost about three times that for computing $e^A$ alone. This is joint work with Awad AlMohy.

10 Oct
2008 Parametric Approximation of Geometric Evolution Equations
John Barrett (Imperial College London)
3.00  Frank Adams Room 1, Alan Turing BuildingAbstract (click to view)Geometric flows, in which hypersurfaces move such that an energy, involving surface and bending terms, decreases appear in many situations in the natural sciences and in geometry. Classic examples are mean curvature, surface diffusion and Willmore flows.
Computational methods to approximate such flows are based on one of three approaches (i) parametric methods, (ii) phase field methods or (iii) level set methods. The first tracks the hypersurface, whilst the other two implicitly capture the hypersurface. A key problem with the first approach, apart from the fact that it is does not naturally deal with changes of topology, is that in many cases the mesh has to be redistributed after every few time steps to avoid coalescence of mesh points.
In this talk we present a variational formulation of the parametric approach, which leads to an unconditionally stable, fully discrete approximation. In addition, the scheme has very good properties with respect to the distribution of mesh points, and if applicable volume conservation. We illustrate this for (anisotropic) mean curvature and (anisotropic) surface diffusion of closed curves in R2. We extend these flows to curve networks in R2. Here the triple junction conditions, that have to hold where three curves meet at a point, are naturally approximated in the discretization of our variational formulation. Our approach naturally generalises to the corresponding flows on curves in Rd, d ≥ 3, as well as to curves on twodimensional manifolds in R3. Finally, we extend these approximations to flows on closed hypersurfaces and to surface clusters in R3. 
24 Oct
2008 Pressure stabilisation for multiphysics problems in fluid mechanics.
Erik Burman (University of Sussex)
3.00  Frank Adams Room 1, Alan Turing BuildingAbstract (click to view)When approximating the Stokes' problem using finite elements it is well known that the velocity and pressure spaces must be chosen in such a way that a certain compatibility condition, the BabuskaBrezzi condition, is satisfied. Velocity/pressure space pairs that do not satisfy the infsup condition may be stabilised using some (weakly) consistent stabilisation operator.
For single phase flows in fixed domain a large number of methods exist in literature and the theory is essentially complete. However more complex applications with multiphysics character are much less well understood. Indeed applications such as for example fluidstructure interaction and stratified twophase flows introduce additional challenges due to moving boundaries across which the model parameters or the model itself changes. The aim of this talk is to show how pressure stabilisation techniques can be used, in combination with Nitschetype domain decomposition, to enhance accuracy or efficiency, for this type of multiphysics problem 
07 Nov
2008 Solving nonlinear semidefinite programming problems by PENNON
Michal Kocvara (University of Birmingham)
3.00  Frank Adams Room 1, Alan Turing BuildingAbstract (click to view)I will present PENNON 1.0  the first release of the code for solving mathematical optimization problems with real and matrix variables, smooth nonlinear (nonconvex) objective and smooth nonlinear equality and inequality constraints. All functions may depend on both types of variables. Additionally, the matrix variables may be subject to spectral constraints, in particular, to positive semidefinite constraints. Matlab and extended AMPL interface allow the user to formulate the problems easily. I will present a few of numerous applications of the code. This is a joint talk with M. Stingl, University of Erlangen.
 14 Nov
2008 A ParameterChoice Method that Exploits Residual Information
Per Christian Hansen (Technical University of Denmark)
3.00  Frank Adams Room 1, Alan Turing BuildingAbstract (click to view)All regularization methods  used for computing stable solutions to inverse problems  rely on the choice of the regularization parameter that balances the solution's smoothness and how well it fits the data. Most algorithms for choosing this regularization parameter are based on the norm of the residual vector. However, the residual vector carries important information about the regularized solution, and therefore it can be advantageous to seek to extract more of the information available there. We present important relations between the residual components and the amount of information that is available in the noisy data, and we show how to use statistical tools and fast Fourier transforms to extract this information efficiently. This approach leads to acomputationally inexpensive parameterchoice rule based on the normalized cumulative periodogram, which is particularly suited for largescale problems.

21 Nov
2008 Adaptive MultiScale Computational Modelling of LargeScale MultiBody Systems with Contact Constraints
Angela Mihai (Cambridge University)
3.00  Frank Adams Room 1, Alan Turing BuildingAbstract (click to view)Constrained partial differential equations (PDEs) modelling static and dynamic multibody contact problems in largescale structural assemblies have a strong impact on various engineering and scientific applications (e.g. automotive and aerospace industries, chemical processing, material science, biotechnology, etc.). In turn, these applications lead to many exciting research problems for which the appropriate mathematical treatment and the development of robust and efficient solution strategies requires the integration of tools from several mathematical subdisciplines (e.g. the theory of optimization and optimal control in a functional analytic setting, the theory of PDEs, numerical analysis, and scientic computing). In this talk, I will present an adaptive multiscale approach for predicting the mechanical behaviour of masonry structures modelled as dynamic frictional multibody contact problems. In this approach, the iterative splitting of the contact problem into normal contact and frictional contact is combined with a semismooth Newton/primaldual activeset procedure to calculate deformations and openings in the model structures. This algorithm is then coupled with a novel adaptive multiscale technique involving a macroscopic (structural) scale and a mesoscopic (individual components) scale to predict appearance of dislocations and stress distribution in largescale masonry structures. Comparisons of the numerical results with data from experimental tests and from practical observations illustrate the predictive capability of the multiscale algorithm. This presentation is based on collaborative research work with Prof. Mark Ainsworth, University of Strathclyde.

5 Dec
2008 Adaptive regularization methods for nonlinear optimization
Coralia Cartis (The University of Edinburgh)
3.00  Frank Adams Room 1, Alan Turing BuildingAbstract (click to view)A new class of methods for nonlinear nonconvex optimization problems will be presented that approximately globally minimizes a quadratic model of the objective regularized by a cubic term, extending to practical largescale problems earlier approaches by Nesterov (2007) and Griewank (1982). Preliminary numerical experiments show our method to perform better than a trustregion implementation, while our convergence and complexity results show it to be at least as reliable to the latter approaches. Extensions to problems with simple constraints will also be presented. Lastly, solving nonlinear leastsquares problems and systems of equations will be discussed, and a lowerorder regularization method will be analysed in this context, namely where a firstorder model of the leastsquares is regularized by a quadratic term. An overestimation property of functions with Lipschitzcontinuous Hessians, and of systems with Lipschitzcontinuous Jacobian underlies and justifies the entire work to be presented. This is joint work with Nick Gould (RAL), Philippe Toint (Namur, Belgium), and for the latter part of my talk, also with Stefania Bellavia and Benedetta Morini (Florence, Italy).

19 Dec
2008 Computing the phifunction of a matrix
Jitse Niesen (University of Leeds)
3.00  Frank Adams Room 1, Alan Turing BuildingAbstract (click to view)The phifunctions are certain functions related to the exponential: the zeroth phifunction is the exponential itself and the first one is defined by $\phi_1(x) = (e^x1) / x$. Phifunctions play an important role in exponential integrators, which are methods for the numerical solution of ordinary differential equations with a stiff linear part and a nonstiff nonlinear part. We will give a brief introduction to exponential integrators and then focus on the computation of the phifunction of a matrix. The central ingredients in the algorithm we present are a Krylovsubspace method for reducing big sparse matrices to a smaller size and a timestepping method akin to the scalingandsquaring method for matrix exponentials. This is joint work with Will Wright (Melbourne).
Semester Two
 13 Feb
2009 A rigorously justified algebraic preconditioner for highcontrast diffusion
Robert Scheichl (University of Bath)
3.00  Frank Adams Room 1, Alan Turing BuildingAbstract (click to view)In this talk I will analyse the robustness of a family of algebraic multigrid (AMG) preconditioners for elliptic PDEs with highcontrast coefficients. Problems with highcontrast coefficients are ubiquitous in porous media flow applications. Consequently, development of efficient solvers for highcontrast heterogeneous media has been a very active area of research. The family of preconditioners considered here exploits the binary character of highcontrast coefficients and is constructed in similar ways as AMG preconditioners by identifying strong and weak couplings in the stiffness matrix. However, for this new family we can rigorously prove robustness and we demonstrate this on a sequence of model problems. Moreover, our numerical experiments show that for sufficiently high contrast the performance of our new preconditioners is almost identical to classical AMG, both in terms of iteration count and CPU time. The theory is based on a singular perturbation analysis.
This work was obtained partly in collaboration with B. Aksoylu (Louisiana State), I.G. Graham (Bath) and H. Klie (ConocoPhillips).  20 Feb
2009 GMRES and oscillatory differential equations
Sheehan Olver (University of Oxford)
3.00  Frank Adams Room 2, Alan Turing BuildingAbstract (click to view)We investigate applying GMRES to differential operators. This method works particularly well when the solutions to the differential equation are oscillatory, and by preconditioning the equation appropriately we can obtain asymptotic decay as the frequency increases. For the case of oscillatory integrals, we will find general conditions for convergence; hence, the method has similar asymptotic properties to the asymptotic expansion, whilst converging for fixed frequencies.
 06 Mar
2009 From sparsity to blocksparsity: solving very large optimization problems with interior point methods
Jacek Gondzio (University of Edinburgh)
3.00  Frank Adams Room 1, Alan Turing BuildingAbstract (click to view)Many reallife economic models involve system dynamics, spatial distribution or uncertainty and lead to largescale optimization problems. Such problems usually have a hidden structure: they are constructed by replication of some small generic block. The linear algebra subproblems which need to be solved by optimization algorithms for such problems involve matrices which are not only sparse, but they additionally display a blockstructure with many smaller blocks sparsely distributed in the large matrix.
I will discuss several advantages of interior point methods and will focus on their implementation which involves the solution of symmetric indefinite system of equations. Next I will discuss how the special blockstructure of very large problems can be exploited in the linear algebra operations. The optimization software OOPS ( ObjectOriented Parallel Solver: http://www.maths.ed.ac.uk/~gondzio/parallel/solver.html ) based on such a technique can efficiently handle very large problems and achieves scalability on a number of different computing platforms.  13 Mar
2009 Parametric approximation of geometric evolution equations and their coupling to
bulk equations
Robert Nürnberg (Imperial College London)
3.00  Frank Adams Room 1, Alan Turing BuildingAbstract (click to view)Geometric flows, in which hypersurfaces move such that an energy, involving surface and bending terms, decreases appear in many situations in the natural sciences and in geometry. Classic examples are mean curvature, surface diffusion and Willmore flows. Computational methods to approximate such flows are based on one of three approaches (i) parametric methods, (ii) phase field methods or (iii) level set methods. The first tracks the hypersurface, whilst the other two implicitly capture the hypersurface. A key problem with the first approach, apart from the fact that it is does not naturally deal with changes of topology, is that in many cases the mesh has to be redistributed after every few time steps to avoid coalescence of mesh points.
In this talk we present a variational formulation of the parametric approach, which leads to an unconditionally stable, fully discrete finite element approximation. In addition, the scheme has very good properties with respect to the distribution of mesh points, and if applicable volume conservation. We illustrate this for (anisotropic) mean curvature and (anisotropic) surface diffusion flows of closed curves in R^2 and closed hypersurfaces in R^3.
Finally, we extend these approximations to the case when the hypersurface motion depends also on an underlying bulk equation, such as the Stefan problem with kinetic undercooling or the MullinsSekerka/HeleShaw flow with surface tension.  27 Mar
2009 Data assimilation in numerical weather prediction: 4DVar and links to other regularisation methods.
Melina Freitag (University of Bath)
3.00  Frank Adams Room 1, Alan Turing BuildingAbstract (click to view)In this talk we will give an introduction to data assimilation techniques as they are used in modern numerical weather prediction. It is wellknown that data assimilation using 4DVar (4D Variation) can be interpreted as some form of Tikhonov regularisation, a very familiar method for solving illposed inverse problems. Such problems appear in a wide range of applications such as geosciences and image restoration, the process of estimating an original image from a given blurred image.
From the latter work it is known that by replacing the $L_2$norm penalty function by an $L_1$norm penalty function the image restoration problems become edgepreserving as they do not penalise the edges of the image. The $L_1$norm penalty regularisation then 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.
This work is supported by GWR (Great Western Research) and the UK MetOffice and joint with N. Nichols (Reading) and C. Budd (Bath).  24 Apr
2009 Smith forms for even and odd matrix polynomials
Christian Mehl (University of Birmingham)
3.00  Frank Adams Room 1, Alan Turing BuildingAbstract (click to view)Even and odd matrix polynomials, i.e., matrix polynomials with coefficient matrices alternating between symmetric and skewsymmetric matrices, frequently arise in applications.
The corresponding eigenvalue problem is usually solved by first linearizing the problem in a structurepreserving way and then by solving the resulting linear even or odd eigenvalue problem.
However, it has been observed that there exist even or odd matrix polynomials that do not allow a structurepreserving linearization.
In this talk, we investigate the possible Smith forms of even and odd matrix polynomials. As a consequence, we will obtain necessary and sufficient condition when a given even or odd matrix
polynomial allows a structurepreserving linearization.
This is joint work with D. Steven Mackey, Niloufer Mackey, and Volker Mehrmann.
 03 Jul 2009 A sprint through the mathematical software landscape
Mike Croucher (University of Manchester)
3.00  Frank Adams Room 1, Alan Turing BuildingAbstract (click to view)The world of mathematical software is becoming increasingly complicated with literally hundreds of different commercial and free products to choose from. In this talk I will discuss the relative merits of the main generalpurpose commercial mathematical applications (and their opensource competitors) available to staff and students at Manchester University.
Further information
For further information please contact the seminar organiser Younes Chahlaoui.