This article describes the current state of the art of high-performance computing systems, and attempts to shed light on near-future developments that might prolong the steady growth in speed of such systems, which has been one of their most remarkable characteristics. We review the different ways devised to speed them up, both with regard to components and their architecture. In addition, we discuss the requirements for software that can take advantage of existing and future architectures.

Retarded functional differential equations (RFDEs) form a wide class of evolution equations which share the property that, at any point, the rate of the solution depends on a discrete or distributed set of values attained by the solution itself in the past. Thus the initial problem for RFDEs is an infinite-dimensional problem, taking its theoretical and numerical analysis beyond the classical schemes developed for differential equations with no functional elements. In particular, numerically solving initial problems for RFDEs is a diffcult task that cannot be founded on the mere adaptation of well-known methods for ordinary, partial or integro-differential equations to the presence of retarded arguments. Indeed, efficient codes for their numerical integration need speciffc approaches designed according to the nature of the equation and the behaviour of the solution.

By defining the numerical method as a suitable approximation of the solution map of the given equation, we present an original and unifying theory for the convergence and accuracy analysis of the approximate solution. Two particular approaches, both inspired by Runge–Kutta methods, are described. Despite being apparently similar, they are intrinsically different. Indeed, in the presence of speciffc types of functionals on the right-hand side, only one of them can have an explicit character, whereas the other gives rise to an overall procedure which is implicit in any case, even for non-stiff problems.

In the panorama of numerical RFDEs, some critical situations have been recently investigated in connection to speciffc classes of equations, such as the accurate location of discontinuity points, the termination and bifurcation of the solutions of neutral equations, with state-dependent delays, the regularization of the equation and the generalization of the solution behind possible termination points, and the treatment of equations stated in the implicit form, which include singularly perturbed problems and delay differential-algebraic equations as well. All these issues are tackled in the last three sections.

In this paper we have not considered the important issue of stability, for which we refer the interested reader to the comprehensive book by Bellen and Zennaro (2003).

Numerical simulation of parametrized differential equations is of crucial importance in the study of real-world phenomena in applied science and engineering. Computational methods for real-time and many-query simulation of such problems often require prohibitively high computational costs to achieve sufficiently accurate numerical solutions. During the last few decades, model order reduction has proved successful in providing low-complexity high-fidelity surrogate models that allow rapid and accurate simulations under parameter variation, thus enabling the numerical simulation of increasingly complex problems. However, many challenges remain to secure the robustness and efficiency needed for the numerical simulation of nonlinear time-dependent problems. The purpose of this article is to survey the state of the art of reduced basis methods for time-dependent problems and draw together recent advances in three main directions. First, we discuss structure-preserving reduced order models designed to retain key physical properties of the continuous problem. Second, we survey localized and adaptive methods based on nonlinear approximations of the solution space. Finally, we consider data-driven techniques based on non-intrusive reduced order models in which an approximation of the map between parameter space and coefficients of the reduced basis is learned. Within each class of methods, we describe different approaches and provide a comparative discussion that lends insights to advantages, disadvantages and potential open questions.

The usual large-scale discretizations are applied to two or three spatial dimensions. The standard methods fail for higher dimensions because the data size increases exponentially with the dimension. In the case of a regular grid with n grid points per direction, a spatial dimension d yields nd grid points. A grid function defined on such a grid is an example of a tensor of order d. Here, suitable tensor formats help, since they try to approximate these huge objects by a much smaller number of parameters, which increases only linearly in d. In this way, data of size nd = 10001000 can also be treated.

This paper introduces the algebraic and analytical aspects of tensor spaces. The main part concerns the numerical representation of tensors and the numerical performance of tensor operations.

Data assimilation addresses the general problem of how to combine model-based predictions with partial and noisy observations of the process in an optimal manner. This survey focuses on sequential data assimilation techniques using probabilistic particle-based algorithms. In addition to surveying recent developments for discrete- and continuous-time data assimilation, both in terms of mathematical foundations and algorithmic implementations, we also provide a unifying framework from the perspective of coupling of measures, and Schrödinger’s boundary value problem for stochastic processes in particular.

This article reviews the application of various notions from the theory of dynamical systems to the analysis of numerical approximation of initial value problems over long-time intervals. Standard error estimates comparing individual trajectories are of no direct use in this context since the error constant typically grows like the exponential of the time interval under consideration.

Instead of comparing trajectories, the effect of discretization on various sets which are invariant under the evolution of the underlying differential equation is studied. Such invariant sets are crucial in determining long-time dynamics. The particular invariant sets which are studied are equilibrium points, together with their unstable manifolds and local phase portraits, periodic solutions, quasi-periodic solutions and strange attractors.

Particular attention is paid to the development of a unified theory and to the development of an existence theory for invariant sets of the underlying differential equation which may be used directly to construct an analogous existence theory (and hence a simple approximation theory) for the numerical method.

In this article we review the image reconstruction algorithms used in tomography. We restrict ourselves to the standard problems in the reconstruction of function from line or plane integrals as they occur in X-ray tomography, nuclear medicine, magnetic resonance imaging, and electron microscopy. Nonstandard situations, such as incomplete data, unknown orientations, local tomography, and discrete tomography are not dealt with. Nor do we treat nonlinear tomographic techniques such as impedance, ultrasound, and near-infrared imaging.

Molecular dynamics is discussed from a mathematical perspective. The recent history of method development is briefly surveyed with an emphasis on the use of geometric integration as a guiding principle. The recovery of statistical mechanical averages from molecular dynamics is then introduced, and the use of backward error analysis as a technique for analysing the accuracy of numerical averages is described. This article gives the first rigorous estimates for the error in statistical averages computed from molecular dynamics simulation based on backward error analysis. It is shown that molecular dynamics introduces an appreciable bias at stepsizes which are below the stability threshold. Simulations performed in such a regime can be corrected by use of a stepsize-dependent reweighting factor. Numerical experiments illustrate the efficacy of this approach. In the final section, several open problems in dynamics-based molecular sampling are considered.

This article starts with a brief introduction to neural networks for those unfamiliar with the basic concepts, together with a very brief overview of mathematical approaches to the subject. This is followed by a more detailed look at three areas of research which are of particular interest to numerical analysts.

The first area is approximation theory. If K is a compact set in ℝn, for some n, then it is proved that a semilinear feedforward network with one hidden layer can uniformly approximate any continuous function in C(K) to any required accuracy. A discussion of known results and open questions on the degree of approximation is included. We also consider the relevance of radial basis functions to neural networks.

The second area considered is that of learning algorithms. A detailed analysis of one popular algorithm (the delta rule) will be given, indicating why one implementation leads to a stable numerical process, whereas an initially attractive variant (essentially a form of steepest descent) does not. Similar considerations apply to the backpropagation algorithm. The effect of filtering and other preprocessing of the input data will also be discussed systematically.

Finally some applications of neural networks to numerical computation are considered.

In the following chapters we will discuss particle methods for the numerical simulation of rarefied gas flows.

We will mainly treat a billiard game, that is, our particles will be hard spheres. But we will also touch upon cases where particles have internal energies due to rotation or vibration, which they exchange in a collision, and we will talk about chemical reactions happening during a collision.

Due to the limited size of this paper, we are only able to mention the principles of these real-gas effects. On the other hand, the general concepts of particle methods to be presented may be used for other kinds of kinetic equations, such as the semiconductor device simulation. We leave this part of the research to subsequent papers.

Let us think about ways to find both eigenvalues and eigenvectors of tridiagonal matrices. An important special case is the computation of singular values and singular vectors of bidiagonal matrices. The discussion is addressed both to specialists in matrix computation and to other scientists whose main interests lie elsewhere. The reason for hoping to communicate with two such diverse sets of readers at the same time is that the content of the survey, though of recent origin, is quite elementary and does not demand familiarity with much beyond triangular factorization and the Gram-Schmidt process for orthogonalizing a set of vectors. For some readers the survey will cover familiar territory but from a novel perspective. The justification for presenting these ideas is that they lead to new variations of current methods that run a lot faster while achieving greater accuracy.

We review level set methods and the related techniques that are common in many PDE-based image models. Many of these techniques involve minimizing the total variation of the solution and admit regularizations on the curvature of its level sets. We examine the scope of these techniques in image science, in particular in image segmentation, interpolation, and decomposition, and introduce some relevant level set techniques that are useful for this class of applications. Many of the standard problems are formulated as variational models. We observe increasing synergistic progression of new tools and ideas between the inverse problem community and the ‘imagers’. We show that image science demands multi-disciplinary knowledge and flexible, but still robust methods. That is why the level set method and total variation methods have become thriving techniques in this field.

Our goal is to survey recently developed techniques in various fields of research that are relevant to diverse objectives in image science. We begin by reviewing some typical PDE-based applications in image processing. In typical PDE methods, images are assumed to be continuous functions sampled on a grid. We will show that these methods all share a common feature, which is the emphasis on processing the level lines of the underlying image. The importance of level lines has been known for some time. See, e.g., Alvarez, Guichard, Morel and Lions (1993). This feature places our slightly general definition of the level set method for image science in context. In Section 2 we describe the building blocks of a typical level set method in the continuum setting. Each important task that we need to do is formulated as the solution to certain PDEs. Then, in Section 3, we briefly describe the finite difference methods developed to construct approximate solutions to these PDEs. Some approaches to interpolation into small subdomains of an image are reviewed in Section 4. In Section 5 we describe the Chan–Vese segmentation algorithm and two new fast implementation methods. Finally, in Section 6, we describe some new techniques developed in the level set community.

In the past decade, problems related to l1/nuclear norm minimization have attracted much attention in the signal processing, machine learning and optimization communities. In this paper, devoted to l1/nuclear norm minimization as ‘optimization beasts’, we give a detailed description of two attractive first-order optimization techniques for solving problems of this type. The first one, aimed primarily at lasso-type problems, comprises fast gradient methods applied to composite minimization formulations. The second approach, aimed at Dantzig-selector-type problems, utilizes saddle-point first-order algorithms and reformulation of the problem of interest as a generalized bilinear saddle-point problem. For both approaches, we give complete and detailed complexity analyses and discuss the application domains.

This article demonstrates how numerical methods for atmospheric models can be validated by showing that they give the theoretically predicted rate of convergence to relevant asymptotic limit solutions. This procedure is necessary because the exact solution of the Navier–Stokes equations cannot be resolved by production models. The limit solutions chosen are those most important for weather and climate prediction. While the best numerical algorithms for this purpose largely reflect current practice, some important limit solutions cannot be captured by existing methods. The use of Lagrangian rather than Eulerian averaging may be required in these cases.

We discuss the modelling framework of port-Hamiltonian descriptor systems and their use in numerical simulation and control. The structure is ideal for automated network-based modelling since it is invariant under power-conserving interconnection, congruence transformations and Galerkin projection. Moreover, stability and passivity properties are easily shown. Condensed forms under orthogonal transformations present easy analysis tools for existence, uniqueness, regularity and numerical methods to check these properties.

After recalling the concepts for general linear and nonlinear descriptor systems, we demonstrate that many difficulties that arise in general descriptor systems can be easily overcome within the port-Hamiltonian framework. The properties of port-Hamiltonian descriptor systems are analysed, and time discretization and numerical linear algebra techniques are discussed. Structure-preserving regularization procedures for descriptor systems are presented to make them suitable for simulation and control. Model reduction techniques that preserve the structure and stabilization and optimal control techniques are discussed.

The properties of port-Hamiltonian descriptor systems and their use in modelling simulation and control methods are illustrated with several examples from different physical domains. The survey concludes with open problems and research topics that deserve further attention.

We describe topological methods for the efficient, rigorous computation of dynamical systems. In particular, we indicate how Conley's Fundamental Decomposition Theorem is naturally related to combinatorial approximations of dynamical systems. Furthermore, we show that computations of Morse decompositions and isolating blocks can be performed efficiently. We conclude with examples indicating how these ideas can be applied to finite- and infinite-dimensional discrete and continuous dynamical systems.

We first survey componentwise and normwise perturbation bounds for the standard least squares (LS) and minimum norm problems. Then some recent estimates of the optimal backward error for an alleged solution to an LS problem are presented. These results are particularly interesting when the algorithm used is not backward stable.

The QR factorization and the singular value decomposition (SVD), developed in the 1960s and early 1970s, remain the basic tools for solving both the LS and the total least squares (TLS) problems. Current algorithms based on Householder or Gram-Schmidt QR factorizations are reviewed. The use of the SVD to determine the numerical rank of a matrix, as well as for computing a sequence of regularized solutions, is then discussed. The solution of the TLS problem in terms of the SVD of the compound matrix $(b\ A)$ is described.

Some recent algorithmic developments are motivated by the need for the efficient implementation of the QR factorization on modern computer architectures. This includes blocked algorithms as well as newer recursive implementations. Other developments come from needs in different application areas. For example, in signal processing rank-revealing orthogonal decompositions need to be frequently updated. We review several classes of such decompositions, which can be more efficiently updated than the SVD.

Two algorithms for the orthogonal bidiagonalization of an arbitrary matrix were given by Golub and Kahan in 1965, one using Householder transformations and the other a Lanczos process. If used to transform the matrix $(b\ A)$ to upper bidiagonal form, this becomes a powerful tool for solving various LS and TLS problems. This bidiagonal decomposition gives a core regular subproblem for the TLS problem. When implemented by the Lanczos process it forms the kernel in the iterative method LSQR. It is also the basis of the partial least squares (PLS) method, which has become a standard tool in statistics.

We present some generalized QR factorizations which can be used to solve different generalized least squares problems. Many applications lead to LS problems where the solution is subject to constraints. This includes linear equality and inequality constraints. Quadratic constraints are used to regularize solutions to discrete ill-posed LS problems. We survey these classes of problems and discuss their solution.

As in all scientific computing, there is a trend that the size and complexity of the problems being solved is steadily growing. Large problems are often sparse or structured. Algorithms for the efficient solution of banded and block-angular LS problems are given, followed by a brief discussion of the general sparse case. Iterative methods are attractive, in particular when matrix-vector multiplication is cheap.

We present the asymptotic transitions from microscopic to macroscopic physics, their computational challenges and the asymptotic-preserving (AP) strategies to compute multiscale physical problems efficiently. Specifically, we will first study the asymptotic transition from quantum to classical mechanics, from classical mechanics to kinetic theory, and then from kinetic theory to hydrodynamics. We then review some representative AP schemes that mimic these asymptotic transitions at the discrete level, and hence can be used crossing scales and, in particular, capture the macroscopic behaviour without resolving the microscopic physical scale numerically.

Kohn–Sham density functional theory (DFT) is the most widely used electronic structure theory. Despite significant progress in the past few decades, the numerical solution of Kohn–Sham DFT problems remains challenging, especially for large-scale systems. In this paper we review the basics as well as state-of-the-art numerical methods, and focus on the unique numerical challenges of DFT.

We provide a ‘user guide’ to the literature of the past twenty years concerning the modelling and approximation of discontinuous solutions to nonlinear hyperbolic systems that admit small-scale dependent shock waves. We cover several classes of problems and solutions: nonclassical undercompressive shocks, hyperbolic systems in nonconservative form, and boundary layer problems. We review the relevant models arising in continuum physics and describe the numerical methods that have been proposed to capture small-scale dependent solutions. In agreement with general well-posedness theory, small-scale dependent solutions are characterized by a kinetic relation, a family of paths, or an admissible boundary set. We provide a review of numerical methods (front-tracking schemes, finite difference schemes, finite volume schemes), which, at the discrete level, reproduce the effect of the physically meaningful dissipation mechanisms of interest in the applications. An essential role is played by the equivalent equation associated with discrete schemes, which is found to be relevant even for solutions containing shock waves.