This paper is concerned with the mathematical structure of the immersed boundary (IB) method, which is intended for the computer simulation of fluid–structure interaction, especially in biological fluid dynamics. The IB formulation of such problems, derived here from the principle of least action, involves both Eulerian and Lagrangian variables, linked by the Dirac delta function. Spatial discretization of the IB equations is based on a fixed Cartesian mesh for the Eulerian variables, and a moving curvilinear mesh for the Lagrangian variables. The two types of variables are linked by interaction equations that involve a smoothed approximation to the Dirac delta function. Eulerian/Lagrangian identities govern the transfer of data from one mesh to the other. Temporal discretization is by a second-order Runge–Kutta method. Current and future research directions are pointed out, and applications of the IB method are briefly discussed.

Large linear systems of saddle point type arise in a wide variety of applications throughout computational science and engineering. Due to their indefiniteness and often poor spectral properties, such linear systems represent a significant challenge for solver developers. In recent years there has been a surge of interest in saddle point problems, and numerous solution techniques have been proposed for this type of system. The aim of this paper is to present and discuss a large selection of solution methods for linear systems in saddle point form, with an emphasis on iterative methods for large and sparse problems.

Since its popularization in the late 1970s, Sequential Quadratic Programming (SQP) has arguably become the most successful method for solving nonlinearly constrained optimization problems. As with most optimization methods, SQP is not a single algorithm, but rather a conceptual method from which numerous specific algorithms have evolved. Backed by a solid theoretical and computational foundation, both commercial and public-domain SQP algorithms have been developed and used to solve a remarkably large set of important practical problems. Recently large-scale versions have been devised and tested with promising results.

The subject of inverse problems in differential equations is of enormous practical importance, and has also generated substantial mathematical and computational innovation. Typically some form of regularization is required to ameliorate ill-posed behaviour. In this article we review the Bayesian approach to regularization, developing a function space viewpoint on the subject. This approach allows for a full characterization of all possible solutions, and their relative probabilities, whilst simultaneously forcing significant modelling issues to be addressed in a clear and precise fashion. Although expensive to implement, this approach is starting to lie within the range of the available computational resources in many application areas. It also allows for the quantification of uncertainty and risk, something which is increasingly demanded by these applications. Furthermore, the approach is conceptually important for the understanding of simpler, computationally expedient approaches to inverse problems.

This paper gives a review of integration algorithms for finite dimensional mechanical systems that are based on discrete variational principles. The variational technique gives a unified treatment of many symplectic schemes, including those of higher order, as well as a natural treatment of the discrete Noether theorem. The approach also allows us to include forces, dissipation and constraints in a natural way. Amongst the many specific schemes treated as examples, the Verlet, SHAKE, RATTLE, Newmark, and the symplectic partitioned Runge–Kutta schemes are presented.

We present a survey of the fundamentals and the applications of sparse grids, with a focus on the solution of partial differential equations (PDEs). The sparse grid approach, introduced in Zenger (1991), is based on a higher-dimensional multiscale basis, which is derived from a one-dimensional multi-scale basis by a tensor product construction. Discretizations on sparse grids involve $O(N \cdot (\log N)^{d-1})$ degrees of freedom only, where $d$ denotes the underlying problem's dimensionality and where $N$ is the number of grid points in one coordinate direction at the boundary. The accuracy obtained with piecewise linear basis functions, for example, is $O(N^{-2} \cdot (\log N)^{d-1})$ with respect to the $L_{2^-}$ and $L_{\infty}$-norm, if the solution has bounded second mixed derivatives. This way, the curse of dimensionality, i.e., the exponential dependence $O(N^d)$ of conventional approaches, is overcome to some extent. For the energy norm, only $O(N)$ degrees of freedom are needed to give an accuracy of $O(N^{-1})$. That is why sparse grids are especially well-suited for problems of very high dimensionality.

The sparse grid approach can be extended to nonsmooth solutions by adaptive refinement methods. Furthermore, it can be generalized from piecewise linear to higher-order polynomials. Also, more sophisticated basis functions like interpolets, prewavelets, or wavelets can be used in a straightforward way.

We describe the basic features of sparse grids and report the results of various numerical experiments for the solution of elliptic PDEs as well as for other selected problems such as numerical quadrature and data mining.

Monte Carlo is one of the most versatile and widely used numerical methods. Its convergence rate, O(N−1/2), is independent of dimension, which shows Monte Carlo to be very robust but also slow. This article presents an introduction to Monte Carlo methods for integration problems, including convergence theory, sampling methods and variance reduction techniques. Accelerated convergence for Monte Carlo quadrature is attained using quasi-random (also called low-discrepancy) sequences, which are a deterministic alternative to random or pseudo-random sequences. The points in a quasi-random sequence are correlated to provide greater uniformity. The resulting quadrature method, called quasi-Monte Carlo, has a convergence rate of approximately O((logN)kN−1). For quasi-Monte Carlo, both theoretical error estimates and practical limitations are presented. Although the emphasis in this article is on integration, Monte Carlo simulation of rarefied gas dynamics is also discussed. In the limit of small mean free path (that is, the fluid dynamic limit), Monte Carlo loses its effectiveness because the collisional distance is much less than the fluid dynamic length scale. Computational examples are presented throughout the text to illustrate the theory. A number of open problems are described.

This article surveys a general approach to error control and adaptive mesh design in Galerkin finite element methods that is based on duality principles as used in optimal control. Most of the existing work on a posteriori error analysis deals with error estimation in global norms like the ‘energy norm’ or the L2 norm, involving usually unknown ‘stability constants’. However, in most applications, the error in a global norm does not provide useful bounds for the errors in the quantities of real physical interest. Further, their sensitivity to local error sources is not properly represented by global stability constants. These deficiencies are overcome by employing duality techniques, as is common in a priori error analysis of finite element methods, and replacing the global stability constants by computationally obtained local sensitivity factors. Combining this with Galerkin orthogonality, a posteriori estimates can be derived directly for the error in the target quantity. In these estimates local residuals of the computed solution are multiplied by weights which measure the dependence of the error on the local residuals. Those, in turn, can be controlled by locally refining or coarsening the computational mesh. The weights are obtained by approximately solving a linear adjoint problem. The resulting a posteriori error estimates provide the basis of a feedback process for successively constructing economical meshes and corresponding error bounds tailored to the particular goal of the computation. This approach, called the ‘dual-weighted-residual method’, is introduced initially within an abstract functional analytic setting, and is then developed in detail for several model situations featuring the characteristic properties of elliptic, parabolic and hyperbolic problems. After having discussed the basic properties of duality-based adaptivity, we demonstrate the potential of this approach by presenting a selection of results obtained for practical test cases. These include problems from viscous fluid flow, chemically reactive flow, elasto-plasticity, radiative transfer, and optimal control. Throughout the paper, open theoretical and practical problems are stated together with references to the relevant literature.

In this paper we consider the construction, analysis, implementation and application of exponential integrators. The focus will be on two types of stiff problems. The first one is characterized by a Jacobian that possesses eigenvalues with large negative real parts. Parabolic partial differential equations and their spatial discretization are typical examples. The second class consists of highly oscillatory problems with purely imaginary eigenvalues of large modulus. Apart from motivating the construction of exponential integrators for various classes of problems, our main intention in this article is to present the mathematics behind these methods. We will derive error bounds that are independent of stiffness or highest frequencies in the system.

Since the implementation of exponential integrators requires the evaluation of the product of a matrix function with a vector, we will briefly discuss some possible approaches as well. The paper concludes with some applications, in which exponential integrators are used.

In this survey we discuss various approximation-theoretic problems that arise in the multilayer feedforward perceptron (MLP) model in neural networks. The MLP model is one of the more popular and practical of the many neural network models. Mathematically it is also one of the simpler models. Nonetheless the mathematics of this model is not well understood, and many of these problems are approximation-theoretic in character. Most of the research we will discuss is of very recent vintage. We will report on what has been done and on various unanswered questions. We will not be presenting practical (algorithmic) methods. We will, however, be exploring the capabilities and limitations of this model.

This is a survey of nonlinear approximation, especially that part of the subject which is important in numerical computation. Nonlinear approximation means that the approximants do not come from linear spaces but rather from nonlinear manifolds. The central question to be studied is what, if any, are the advantages of nonlinear approximation over the simpler, more established, linear methods. This question is answered by studying the rate of approximation which is the decrease in error versus the number of parameters in the approximant. The number of parameters usually correlates well with computational effort. It is shown that in many settings the rate of nonlinear approximation can be characterized by certain smoothness conditions which are significantly weaker than required in the linear theory. Emphasis in the survey will be placed on approximation by piecewise polynomials and wavelets as well as their numerical implementation. Results on highly nonlinear methods such as optimal basis selection and greedy algorithms (adaptive pursuit) are also given. Applications to image processing, statistical estimation, regularity for PDEs, and adaptive algorithms are discussed.

Finite element exterior calculus is an approach to the design and understanding of finite element discretizations for a wide variety of systems of partial differential equations. This approach brings to bear tools from differential geometry, algebraic topology, and homological algebra to develop discretizations which are compatible with the geometric, topological, and algebraic structures which underlie well-posedness of the PDE problem being solved. In the finite element exterior calculus, many finite element spaces are revealed as spaces of piecewise polynomial differential forms. These connect to each other in discrete subcomplexes of elliptic differential complexes, and are also related to the continuous elliptic complex through projections which commute with the complex differential. Applications are made to the finite element discretization of a variety of problems, including the Hodge Laplacian, Maxwell’s equations, the equations of elasticity, and elliptic eigenvalue problems, and also to preconditioners.

We introduce a new version of the Fast Multipole Method for the evaluation of potential fields in three dimensions. It is based on a new diagonal form for translation operators and yields high accuracy at a reasonable cost.

Radial basis function methods are modern ways to approximate multivariate functions, especially in the absence of grid data. They have been known, tested and analysed for several years now and many positive properties have been identified. This paper gives a selective but up-to-date survey of several recent developments that explains their usefulness from the theoretical point of view and contributes useful new classes of radial basis function. We consider particularly the new results on convergence rates of interpolation with radial basis functions, as well as some of the various achievements on approximation on spheres, and the efficient numerical computation of interpolants for very large sets of data. Several examples of useful applications are stated at the end of the paper.

This article discusses finite element Galerkin schemes for a number of linear model problems in electromagnetism. The finite element schemes are introduced as discrete differential forms, matching the coordinate-independent statement of Maxwell's equations in the calculus of differential forms. The asymptotic convergence of discrete solutions is investigated theoretically. As discrete differential forms represent a genuine generalization of conventional Lagrangian finite elements, the analysis is based upon a judicious adaptation of established techniques in the theory of finite elements. Risks and difficulties haunting finite element schemes that do not fit the framework of discrete differential forms are highlighted.

I thought that instead of the great number of precepts of which logic is composed, I would have enough with the four following ones, provided that I made a firm and unalterable resolution not to violate them even in a single instance. The first rule was never to accept anything as true unless I recognized it to be certainly and evidently such …. The second was to divide each of the difficulties which I encountered into as many parts as possible, and as might be required for an easier solution. (Descartes)

We survey splitting methods for the numerical integration of ordinary differential equations (ODEs). Splitting methods arise when a vector field can be split into a sum of two or more parts that are each simpler to integrate than the original (in a sense to be made precise). One of the main applications of splitting methods is in geometric integration, that is, the integration of vector fields that possess a certain geometric property (e.g., being Hamiltonian, or divergence-free, or possessing a symmetry or first integral) that one wants to preserve. We first survey the classification of geometric properties of dynamical systems, before considering the theory and applications of splitting in each case. Once a splitting is constructed, the pieces are composed to form the integrator; we discuss the theory of such ‘composition methods’ and summarize the best currently known methods. Finally, we survey applications from celestial mechanics, quantum mechanics, accelerator physics, molecular dynamics, and fluid dynamics, and examples from dynamical systems, biology and reaction–diffusion systems.

Many differential equations of practical interest evolve on Lie groups or on manifolds acted upon by Lie groups. The retention of Lie-group structure under discretization is often vital in the recovery of qualitatively correct geometry and dynamics and in the minimization of numerical error. Having introduced requisite elements of differential geometry, this paper surveys the novel theory of numerical integrators that respect Lie-group structure, highlighting theory, algorithmic issues and a number of applications.

Many optimal decision problems in scientific, engineering, and public sector applications involve both discrete decisions and nonlinear system dynamics that affect the quality of the final design or plan. These decision problems lead to mixed-integer nonlinear programming (MINLP) problems that combine the combinatorial difficulty of optimizing over discrete variable sets with the challenges of handling nonlinear functions. We review models and applications of MINLP, and survey the state of the art in methods for solving this challenging class of problems.

Most solution methods for MINLP apply some form of tree search. We distinguish two broad classes of methods: single-tree and multitree methods. We discuss these two classes of methods first in the case where the underlying problem functions are convex. Classical single-tree methods include nonlinear branch-and-bound and branch-and-cut methods, while classical multitree methods include outer approximation and Benders decomposition. The most efficient class of methods for convex MINLP are hybrid methods that combine the strengths of both classes of classical techniques.

Non-convex MINLPs pose additional challenges, because they contain non-convex functions in the objective function or the constraints; hence even when the integer variables are relaxed to be continuous, the feasible region is generally non-convex, resulting in many local minima. We discuss a range of approaches for tackling this challenging class of problems, including piecewise linear approximations, generic strategies for obtaining convex relaxations for non-convex functions, spatial branch-and-bound methods, and a small sample of techniques that exploit particular types of non-convex structures to obtain improved convex relaxations.

We finish our survey with a brief discussion of three important aspects of MINLP. First, we review heuristic techniques that can obtain good feasible solution in situations where the search-tree has grown too large or we require real-time solutions. Second, we describe an emerging area of mixed-integer optimal control that adds systems of ordinary differential equations to MINLP. Third, we survey the state of the art in software for MINLP.

Monte Carlo methods are a very general and useful approach for the estimation of expectations arising from stochastic simulation. However, they can be computationally expensive, particularly when the cost of generating individual stochastic samples is very high, as in the case of stochastic PDEs. Multilevel Monte Carlo is a recently developed approach which greatly reduces the computational cost by performing most simulations with low accuracy at a correspondingly low cost, with relatively few simulations being performed at high accuracy and a high cost.

In this article, we review the ideas behind the multilevel Monte Carlo method, and various recent generalizations and extensions, and discuss a number of applications which illustrate the flexibility and generality of the approach and the challenges in developing more efficient implementations with a faster rate of convergence of the multilevel correction variance.

The subject of geometric numerical integration deals with numerical integrators that preserve geometric properties of the flow of a differential equation, and it explains how structure preservation leads to improved long-time behaviour. This article illustrates concepts and results of geometric numerical integration on the important example of the Störmer–Verlet method. It thus presents a cross-section of the recent monograph by the authors, enriched by some additional material.

After an introduction to the Newton–Störmer–Verlet–leapfrog method and its various interpretations, there follows a discussion of geometric properties: reversibility, symplecticity, volume preservation, and conservation of first integrals. The extension to Hamiltonian systems on manifolds is also described. The theoretical foundation relies on a backward error analysis, which translates the geometric properties of the method into the structure of a modified differential equation, whose flow is nearly identical to the numerical method. Combined with results from perturbation theory, this explains the excellent long-time behaviour of the method: long-time energy conservation, linear error growth and preservation of invariant tori in near-integrable systems, a discrete virial theorem, and preservation of adiabatic invariants.