In the sciences, situations where dissipation is not significant may invariably be modelled by Hamiltonian systems of ordinary, or partial, differential equations. Symplectic integrators are numerical methods specifically aimed at advancing in time the solution of Hamiltonian systems. Roughly speaking, ‘symplecticness’ is a characteristic property possessed by the solutions of Hamiltonian problems. A numerical method is called symplectic if, when applied to Hamiltonian problems, it generates numerical solutions which inherit the property of symplecticness.

In recent years, reduced-order modelling techniques based on Krylov-subspace iterations, especially the Lanczos algorithm and the Arnoldi process, have become popular tools for tackling the large-scale time-invariant linear dynamical systems that arise in the simulation of electronic circuits. This paper reviews the main ideas of reduced-order modelling techniques based on Krylov subspaces and describes some applications of reduced-order modelling in circuit simulation.

Regularization methods are a key tool in the solution of inverse problems. They are used to introduce prior knowledge and allow a robust approximation of ill-posed (pseudo-) inverses. In the last two decades interest has shifted from linear to nonlinear regularization methods, even for linear inverse problems. The aim of this paper is to provide a reasonably comprehensive overview of this shift towards modern nonlinear regularization methods, including their analysis, applications and issues for future research.

In particular we will discuss variational methods and techniques derived from them, since they have attracted much recent interest and link to other fields, such as image processing and compressed sensing. We further point to developments related to statistical inverse problems, multiscale decompositions and learning theory.

In many optimization problems arising from scientific, engineering and artificial intelligence applications, objective and constraint functions are available only as the output of a black-box or simulation oracle that does not provide derivative information. Such settings necessitate the use of methods for derivative-free, or zeroth-order, optimization. We provide a review and perspectives on developments in these methods, with an emphasis on highlighting recent developments and on unifying treatment of such problems in the non-linear optimization and machine learning literature. We categorize methods based on assumed properties of the black-box functions, as well as features of the methods. We first overview the primary setting of deterministic methods applied to unconstrained, non-convex optimization problems where the objective function is defined by a deterministic black-box oracle. We then discuss developments in randomized methods, methods that assume some additional structure about the objective (including convexity, separability and general non-smooth compositions), methods for problems where the output of the black-box oracle is stochastic, and methods for handling different types of constraints.

This review paper collects several results that form part of the theoretical foundation of isogeometric methods. We analyse variational techniques for the numerical resolution of PDEs based on splines or NURBS and we provide optimal approximation and error estimates in several cases of interest. The theory presented also includes estimates for T-splines, which are an extension of splines allowing for local refinement. In particular, we focus our attention on elliptic and saddle point problems, and we define spline edge and face elements. Our theoretical results are demonstrated by a rich set of numerical examples. Finally, we discuss implementation and efficiency together with preconditioning issues for the final linear system.

Finite differences provided the first numerical approach that permitted large-scale simulations in many applications areas, such as geophysical fluid dynamics. As accuracy and integration time requirements gradually increased, the focus shifted from finite differences to a variety of different spectral methods. During the last few years, radial basis functions, in particular in their ‘local’ RBF-FD form, have taken the major step from being mostly a curiosity approach for small-scale PDE ‘toy problems’ to becoming a major contender also for very large simulations on advanced distributed memory computer systems. Being entirely mesh-free, RBF-FD discretizations are also particularly easy to implement, even when local refinements are needed. This article gives some background to this development, and highlights some recent results.

Numerical modelling of transoceanic tsunami propagation, together with the detailed modelling of inundation of small-scale coastal regions, poses a number of algorithmic challenges. The depth-averaged shallow water equations can be used to reduce this to a time-dependent problem in two space dimensions, but even so it is crucial to use adaptive mesh refinement in order to efficiently handle the vast differences in spatial scales. This must be done in a ‘wellbalanced’ manner that accurately captures very small perturbations to the steady state of the ocean at rest. Inundation can be modelled by allowing cells to dynamically change from dry to wet, but this must also be done carefully near refinement boundaries. We discuss these issues in the context of Riemann-solver-based finite volume methods for tsunami modelling. Several examples are presented using the GeoClaw software, and sample codes are available to accompany the paper. The techniques discussed also apply to a variety of other geophysical flows.

A classical mathematical proof is constructed using pencil and paper. However, there are many ways in which computers may be used in a mathematical proof. But ‘proof by computer’, or even the use of computers in the course of a proof, is not so readily accepted (the December 2008 issue of the Notices of the American Mathematical Society is devoted to formal proofs by computer).

In the following we introduce verification methods and discuss how they can assist in achieving a mathematically rigorous result. In particular we emphasize how floating-point arithmetic is used.

Recent advances in the field of iterative methods for solving large linear systems are reviewed. The main focus is on developments in the area of conjugate gradient-type algorithms and Krylov subspace methods for nonHermitian matrices.

This survey covers the state of the art of techniques for solving general-purpose constrained global optimization problems and continuous constraint satisfaction problems, with emphasis on complete techniques that provably find all solutions (if there are finitely many). The core of the material is presented in sufficient detail that the survey may serve as a text for teaching constrained global optimization.

After giving motivations for and important examples of applications of global optimization, a precise problem definition is given, and a general form of the traditional first-order necessary conditions for a solution. Then more than a dozen software packages for complete global search are described.

A quick review of incomplete methods for bound-constrained problems and recipes for their use in the constrained case follows; an explicit example is discussed, introducing the main techniques used within branch and bound techniques. Sections on interval arithmetic, constrained propagation and local optimization are followed by a discussion of how to avoid the cluster problem. Then a discussion of important problem transformations follows, in particular of linear, convex, and semilinear (= mixed integer linear) relaxations that are important for handling larger problems.

Next, reliability issues – centring on rounding error handling and testing methodologies – are discussed, and the COCONUT framework for the integration of the different techniques is introduced. A list of challenges facing the field in the near future concludes the survey.

Kernels are valuable tools in various fields of numerical analysis, including approximation, interpolation, meshless methods for solving partial differential equations, neural networks, and machine learning. This contribution explains why and how kernels are applied in these disciplines. It uncovers the links between them, in so far as they are related to kernel techniques. It addresses non-expert readers and focuses on practical guidelines for using kernels in applications.

In this article we survey r-adaptive (or moving grid) methods for solving time-dependent partial differential equations (PDEs). Although these methods have received much less attention than their h- and p-adaptive counterparts, particularly within the finite element community, we review the substantial progress that has been made in developing more robust and reliable algorithms and in understanding the basic principles behind these methods, and we give some numerical examples illustrative of the wide classes of problems for which these methods are suitable alternatives to the traditional ones.

More specifically, we first examine the basic geometric properties of moving meshes in both one and higher spatial dimensions, and discuss the discretization process for PDEs on such moving meshes (both structured and unstructured). In particular, we consider the issues of mesh regularity, equidistribution, alignment, and associated variational methods. An overview is given of the general interpolation error analysis for a function or a truncation error on such an adaptive mesh. Guided by these principles, we show how to design effective moving mesh strategies. We then examine in more detail how these strategies can be implemented in practice. The first class of methods which we consider are based upon controlling mesh density and hence are called position-based methods. These make use of a so-called moving mesh PDE (MMPDE) approach and variational methods, as well as optimal transport methods. This is followed by an analysis of methods which have a more Lagrange-like interpretation, and due to this focus are called velocity-based methods. These include the moving finite element method (MFE), the geometric conservation law (GCL) methods, and the deformation map method. Finally, we present a number of specific types of examples for which the use of a moving mesh method is particularly effective in applications. These include scale-invariant problems, blow-up problems, problems with moving fronts and problems in meteorology. We conclude that, whilst r-adaptive methods are still in their relatively early stages of development, with many outstanding questions remaining, they have enormous potential and indeed can produce an optimal form of adaptivity for many problems.

A few months ago, while preparing a lecture to an audience that included engineers and numerical analysts, I asked myself the question: from the point of view of a user of nonlinear optimization routines, how interesting and practical is the body of theoretical analysis developed in this field? To make the question a bit more precise, I decided to select the best optimization methods known to date – those methods that deserve to be in a subroutine library – and for each method ask: what do we know about the behaviour of this method, as implemented in practice? To make my task more tractable, I decided to consider only algorithms for unconstrained optimization.

Partial differential equations (PDEs) with random input data, such as random loadings and coefficients, are reformulated as parametric, deterministic PDEs on parameter spaces of high, possibly infinite dimension. Tensorized operator equations for spatial and temporal k-point correlation functions of their random solutions are derived. Parametric, deterministic PDEs for the laws of the random solutions are derived. Representations of the random solutions' laws on infinite-dimensional parameter spaces in terms of ‘generalized polynomial chaos’ (GPC) series are established. Recent results on the regularity of solutions of these parametric PDEs are presented. Convergence rates of best N-term approximations, for adaptive stochastic Galerkin and collocation discretizations of the parametric, deterministic PDEs, are established. Sparse tensor products of hierarchical (multi-level) discretizations in physical space (and time), and GPC expansions in parameter space, are shown to converge at rates which are independent of the dimension of the parameter space. A convergence analysis of multi-level Monte Carlo (MLMC) discretizations of PDEs with random coefficients is presented. Sufficient conditions on the random inputs for superiority of sparse tensor discretizations over MLMC discretizations are established for linear elliptic, parabolic and hyperbolic PDEs with random coefficients.

Optimization problems involving eigenvalues arise in many different mathematical disciplines. This article is divided into two parts. Part I gives a historical account of the development of the field. We discuss various applications that have been especially influential, from structural analysis to combinatorial optimization, and we survey algorithmic developments, including the recent advance of interior-point methods for a specific problem class: semidefinite programming. In Part II we primarily address optimization of convex functions of eigenvalues of symmetric matrices subject to linear constraints. We derive a fairly complete mathematical theory, some of it classical and some of it new. Using the elegant language of conjugate duality theory, we highlight the parallels between the analysis of invariant matrix norms and weakly invariant convex matrix functions. We then restrict our attention further to linear and semidefinite programming, emphasizing the parallel duality theory and comparing primal-dual interior-point methods for the two problem classes. The final section presents some apparently new variational results about eigenvalues of nonsymmetric matrices, unifying known characterizations of the spectral abscissa (related to Lyapunov theory) and the spectral radius (as an infimum of matrix norms).

This paper aims to give an overview and summary of numerical methods for the solution of stochastic differential equations. It covers discrete time strong and weak approximation methods that are suitable for different applications. A range of approaches and results is discussed within a unified framework. On the one hand, these methods can be interpreted as generalizing the well-developed theory on numerical analysis for deterministic ordinary differential equations. On the other hand they highlight the specific stochastic nature of the equations. In some cases these methods lead to completely new and challenging problems.

In this article we describe recent progress on the design, analysis and implementation of hybrid numerical-asymptotic boundary integral methods for boundary value problems for the Helmholtz equation that model time harmonic acoustic wave scattering in domains exterior to impenetrable obstacles. These hybrid methods combine conventional piecewise polynomial approximations with high-frequency asymptotics to build basis functions suitable for representing the oscillatory solutions. They have the potential to solve scattering problems accurately in a computation time that is (almost) independent of frequency and this has been realized for many model problems. The design and analysis of this class of methods requires new results on the analysis and numerical analysis of highly oscillatory boundary integral operators and on the high-frequency asymptotics of scattering problems. The implementation requires the development of appropriate quadrature rules for highly oscillatory integrals. This article contains a historical account of the development of this currently very active field, a detailed account of recent progress and, in addition, a number of original research results on the design, analysis and implementation of these methods.

Parametrized families of PDEs arise in various contexts such as inverse problems, control and optimization, risk assessment, and uncertainty quantification. In most of these applications, the number of parameters is large or perhaps even infinite. Thus, the development of numerical methods for these parametric problems is faced with the possible curse of dimensionality. This article is directed at (i) identifying and understanding which properties of parametric equations allow one to avoid this curse and (ii) developing and analysing effective numerical methods which fully exploit these properties and, in turn, are immune to the growth in dimensionality.

Part I of this article studies the smoothness and approximability of the solution map, that is, the map $a\mapsto u(a)$, where $a$ is the parameter value and $u(a)$ is the corresponding solution to the PDE. It is shown that for many relevant parametric PDEs, the parametric smoothness of this map is typically holomorphic and also highly anisotropic, in that the relevant parameters are of widely varying importance in describing the solution. These two properties are then exploited to establish convergence rates of $n$-term approximations to the solution map, for which each term is separable in the parametric and physical variables. These results reveal that, at least on a theoretical level, the solution map can be well approximated by discretizations of moderate complexity, thereby showing how the curse of dimensionality is broken. This theoretical analysis is carried out through concepts of approximation theory such as best $n$-term approximation, sparsity, and $n$-widths. These notions determine a priori the best possible performance of numerical methods and thus serve as a benchmark for concrete algorithms.

Part II of this article turns to the development of numerical algorithms based on the theoretically established sparse separable approximations. The numerical methods studied fall into two general categories. The first uses polynomial expansions in terms of the parameters to approximate the solution map. The second one searches for suitable low-dimensional spaces for simultaneously approximating all members of the parametric family. The numerical implementation of these approaches is carried out through adaptive and greedy algorithms. An a priori analysis of the performance of these algorithms establishes how well they meet the theoretical benchmarks.

Mathematical and numerical modelling of the cardiovascular system is a research topic that has attracted remarkable interest from the mathematical community because of its intrinsic mathematical difficulty and the increasing impact of cardiovascular diseases worldwide. In this review article we will address the two principal components of the cardiovascular system: arterial circulation and heart function. We will systematically describe all aspects of the problem, ranging from data imaging acquisition, stating the basic physical principles, analysing the associated mathematical models that comprise PDE and ODE systems, proposing sound and efficient numerical methods for their approximation, and simulating both benchmark problems and clinically inspired problems. Mathematical modelling itself imposes tremendous challenges, due to the amazing complexity of the cardiocirculatory system, the multiscale nature of the physiological processes involved, and the need to devise computational methods that are stable, reliable and efficient. Critical issues involve filtering the data, identifying the parameters of mathematical models, devising optimal treatments and accounting for uncertainties. For this reason, we will devote the last part of the paper to control and inverse problems, including parameter estimation, uncertainty quantification and the development of reduced-order models that are of paramount importance when solving problems with high complexity, which would otherwise be out of reach.

We consider variationally consistent discretization schemes for mechanical contact problems. Most of the results can also be applied to other variational inequalities, such as those for phase transition problems in porous media, for plasticity or for option pricing applications from finance. The starting point is to weakly incorporate the constraint into the setting and to reformulate the inequality in the displacement in terms of a saddle-point problem. Here, the Lagrange multiplier represents the surface forces, and the constraints are restricted to the boundary of the simulation domain. Having a uniform inf-sup bound, one can then establish optimal low-order a priori convergence rates for the discretization error in the primal and dual variables. In addition to the abstract framework of linear saddle-point theory, complementarity terms have to be taken into account. The resulting inequality system is solved by rewriting it equivalently by means of the non-linear complementarity function as a system of equations. Although it is not differentiable in the classical sense, semi-smooth Newton methods, yielding super-linear convergence rates, can be applied and easily implemented in terms of a primal–dual active set strategy. Quite often the solution of contact problems has a low regularity, and the efficiency of the approach can be improved by using adaptive refinement techniques. Different standard types, such as residual- and equilibrated-based a posteriori error estimators, can be designed based on the interpretation of the dual variable as Neumann boundary condition. For the fully dynamic setting it is of interest to apply energy-preserving time-integration schemes. However, the differential algebraic character of the system can result in high oscillations if standard methods are applied. A possible remedy is to modify the fully discretized system by a local redistribution of the mass. Numerical results in two and three dimensions illustrate the wide range of possible applications and show the performance of the space discretization scheme, non-linear solver, adaptive refinement process and time integration.