In this paper we combine the theory of reproducing kernel Hilbert spaces with the field of collocation methods to solve boundary value problems with special emphasis on reproducing property of kernels. From the reproducing property of kernels we proposed a new efficient algorithm to obtain the cardinal functions of a reproducing kernel Hilbert space which can be apply conveniently for multidimensional domains. The differentiation matrices are constructed and also we drive pointwise error estimate of applying them. In addition we prove the nonsingularity of collocation matrix. The proposed method is truly meshless and can be applied conveniently and accurately for high order and also multidimensional problems. Numerical results are presented for the several problems such as second and fifth order two point boundary value problems, one and two dimensional unsteady Burgers equations and a parabolic partial differential equation in three dimensions. Also we compare the numerical results with those reported in the literature to show the high accuracy and efficiency of the proposed method

Dense kernel matrices $\Theta \in \mathbb{R}^{N \times N}$ obtained from point evaluations of a covariance function $G$ at locations $\{ x_{i} \}_{1 \leq i \leq N} \subset \mathbb{R}^{d}$ arise in statistics, machine learning, and numerical analysis. For covariance functions that are Green's functions of elliptic boundary value problems and homogeneously-distributed sampling points, we show how to identify a subset $S \subset \{ 1 , \dots , N \}^2$, with $\# S = O ( N \log (N) \log^{d} ( N /\epsilon ) )$, such that the zero fill-in incomplete Cholesky factorisation of the sparse matrix $\Theta_{ij} 1_{( i, j ) \in S}$ is an $\epsilon$-approximation of $\Theta$. This factorisation can provably be obtained in complexity $O ( N \log( N ) \log^{d}( N /\epsilon) )$ in space and $O ( N \log^{2}( N ) \log^{2d}( N /\epsilon) )$ in time, improving upon the state of the art for general elliptic operators; we further present numerical evidence that $d$ can be taken to be the intrinsic dimension of the data set rather than that of the ambient space. The algorithm only needs to know the spatial configuration of the $x_{i}$ and does not require an analytic representation of $G$. Furthermore, this factorization straightforwardly provides an approximate sparse PCA with optimal rate of convergence in the operator norm. Hence, by using only subsampling and the incomplete Cholesky factorization, we obtain, at nearly linear complexity, the compression, inversion and approximate PCA of a large class of covariance matrices. By inverting the order of the Cholesky factorization we also obtain a solver for elliptic PDE with complexity $O ( N \log^{d}( N /\epsilon) )$ in space and $O ( N \log^{2d}( N /\epsilon) )$ in time, improving upon the state of the art for general elliptic operators.

Lane changing is one of the most common maneuvers on motorways. Although, macroscopic traffic models are well known for their suitability to describe fast moving crowded traffic, most of these models are generally developed in one dimensional framework, henceforth lane changing behavior is somehow neglected. In this paper, we propose a macroscopic model, which accounts for lane-changing behavior on motorway, based on a two-dimensional extension of the Aw and Rascle [Aw and Rascle, SIAM J.Appl.Math., 2000] and Zhang [Zhang, Transport.Res.B-Meth., 2002] macroscopic model for traffic flow. Under conditions, when lane changing maneuvers are no longer possible, the model "relaxes" to the one-dimensional Aw-Rascle-Zhang model. Following the same approach as in [Aw, Klar, Materne and Rascle, SIAM J.Appl.Math., 2002], we derive the two-dimensional macroscopic model through scaling of time discretization of a microscopic follow-the-leader model with driving direction. We provide a detailed analysis of the space-time discretization of the proposed macroscopic as well as an approximation of the solution to the associated Riemann problem. Furthermore, we illustrate some features of the proposed model through some numerical experiments.

Wavelet (Besov) priors are a promising way of reconstructing indirectly measured fields in a regularized manner. We demonstrate how wavelets can be used as a localized basis for reconstructing permeability fields with sharp interfaces from noisy pointwise pressure field measurements in the context of the elliptic inverse problem. For this we derive the adjoint method of minimizing the Besov-norm-regularized misfit functional (this corresponds to determining the maximum a posteriori point in the Bayesian point of view) in the Haar wavelet setting. As it turns out, choosing a wavelet--based prior allows for accelerated optimization compared to established trigonometrically--based priors.

Optimization over low rank matrices has broad applications in machine learning. For large scale problems, an attractive heuristic is to factorize the low rank matrix to a product of two much smaller matrices. In this paper, we study the nonconvex problem $\min_{U\in\mathcal{R}^{n\times r}} g(U)=f(UU^T)$ under the assumptions that $f(X)$ is restricted $\mu$-strongly convex and $L$-smooth on the set $\{X:X\succeq 0,rank(X)\leq r\}$. We propose an accelerated gradient method with alternating constraint that operates directly on the $U$ factors and show that the method has local linear convergence rate with the optimal dependence on the condition number of $\sqrt{L/\mu}$. Globally, our method converges to the critical point with zero gradient from any initializer. Our method also applies to the problem with the asymmetric factorization of $X=\widetilde U\widetilde V^T$ and the same convergence result can be obtained. Extensive experimental results verify the advantage of our method.

This paper investigates gradient recovery schemes for data defined on discretized manifolds. The proposed method, parametric polynomial preserving recovery (PPPR), does not require the tangent spaces of the exact manifolds, and they have been assumed for some significant gradient recovery methods in the literature. Another advantage of PPPR is that superconvergence is guaranteed without the symmetric condition which has been asked in the existing techniques. There is also numerical evidence that the superconvergence by PPPR is high curvature stable, which distinguishes itself from the others. As an application, we show its capability of constructing an asymptotically exact \textit{a posteriori} error estimator. Several numerical examples on two-dimensional surfaces are presented to support the theoretical results and comparisons with existing methods are documented.

We study the exponent of the exponential rate of convergence in terms of the number of degrees of freedom for various non-standard {$p$-version} finite element spaces employing reduced cardinality basis. More specifically, we show that serendipity finite element methods and discontinuous Galerkin finite element methods with total degree $\mathcal{P}_p$ basis have a faster exponential convergence with respect to the number of degrees of freedom than their counterparts employing the tensor product $\mathcal{Q}_p$ basis for quadrilateral/hexahedral elements, for piecewise analytic problems under $p$-refinement. The above results are proven by using a new $p$-optimal error bound for the $L^2$-orthogonal projection onto the total degree $\mathcal{P}_p$ basis, and for the $H^1$-projection onto the serendipity finite element space over tensor product elements with dimension $d\geq2$. These new $p$-optimal error bounds lead to a larger exponent of the exponential rate of convergence with respect to the number of degrees of freedom. Moreover, these results show that part of the basis functions in $\mathcal{Q}_p$ basis {plays} no roles in achieving the $hp$-optimal error bound in the Sobolev space. The sharpness of theoretical results is also verified by a series of numerical examples.

Algebraic multigrid (AMG) is often an effective solver for symmetric positive definite (SPD) linear systems resulting from the discretization of general elliptic PDEs, or the spatial discretization of parabolic PDEs. However, convergence theory and most variations of AMG rely on $A$ being SPD. Hyperbolic PDEs, which arise often in large-scale scientific simulations, remain a challenge for AMG, as well as other fast linear solvers, in part because the resulting linear systems are often highly nonsymmetric. Here, a novel convergence framework is developed for nonsymmetric, reduction-based AMG, and sufficient conditions derived for $\ell^2$-convergence of error and residual. In particular, classical multigrid approximation properties are connected with reduction-based measures to develop a robust framework for nonsymmetric, reduction-based AMG. Matrices with block-triangular structure are then recognized as being amenable to reduction-type algorithms, and a reduction-based AMG method is developed for upwind discretizations of hyperbolic PDEs, based on the concept of a Neumann approximation to ideal restriction ($n$AIR). $n$AIR can be seen as a variation of local AIR ($\ell$AIR) introduced in previous work, specifically targeting matrices with triangular structure. Although less versatile than $\ell$AIR, setup times for $n$AIR can be substantially faster for problems with high connectivity. $n$AIR is shown to be an effective and scalable solver of steady state transport for discontinuous, upwind discretizations, with unstructured meshes, and up to 6th-order finite elements, offering a significant improvement over existing AMG methods. $n$AIR is also shown to be effective on several classes of `nearly triangular' matrices, resulting from curvilinear finite elements and artificial diffusion.

This paper presents a weakly intrusive strategy for computing a low-rank approximation of the solution of a system of nonlinear parameter-dependent equations. The proposed strategy relies on a Newton-like iterative solver which only requires evaluations of the residual of the parameter-dependent equation and of a preconditioner (such as the differential of the residual) for instances of the parameters independently. The algorithm provides an approximation of the set of solutions associated with a possibly large number of instances of the parameters, with a computational complexity which can be orders of magnitude lower than when using the same Newton-like solver for all instances of the parameters. The reduction of complexity requires efficient strategies for obtaining low-rank approximations of the residual, of the preconditioner, and of the increment at each iteration of the algorithm. For the approximation of the residual and the preconditioner, weakly intrusive variants of the empirical interpolation method are introduced, which require evaluations of entries of the residual and the preconditioner. Then, an approximation of the increment is obtained by using a greedy algorithm for low-rank approximation, and a low-rank approximation of the iterate is finally obtained by using a truncated singular value decomposition. When the preconditioner is the differential of the residual, the proposed algorithm is interpreted as an inexact Newton solver for which a detailed convergence analysis is provided. Numerical examples illustrate the efficiency of the method.

As a counterpoint to classical stochastic particle methods for diffusion, we develop a deterministic particle method for linear and nonlinear diffusion. At first glance, deterministic particle methods are incompatible with diffusive partial differential equations since initial data given by sums of Dirac masses would be smoothed instantaneously: particles do not remain particles. Inspired by classical vortex blob methods, we introduce a nonlocal regularization of our velocity field that ensures particles do remain particles, and we apply this to develop a numerical blob method for a range of diffusive partial differential equations of Wasserstein gradient flow type, including the heat equation, the porous medium equation, the Fokker-Planck equation, the Keller-Segel equation, and its variants. Our choice of regularization is guided by the Wasserstein gradient flow structure, and the corresponding energy has a novel form, combining aspects of the well-known interaction and potential energies. In the presence of a confining drift or interaction potential, we prove that minimizers of the regularized energy exist and, as the regularization is removed, converge to the minimizers of the unregularized energy. We then restrict our attention to nonlinear diffusion of porous medium type with at least quadratic exponent. Under sufficient regularity assumptions, we prove that gradient flows of the regularized energies converge to solutions of the porous medium equation. As a corollary, we obtain convergence of our numerical blob method, again under sufficient regularity assumptions. We conclude by considering a range of numerical examples to demonstrate our method's rate of convergence to exact solutions and to illustrate key qualitative properties preserved by the method, including asymptotic behavior of the Fokker-Planck equation and critical mass of the two-dimensional Keller-Segel equation.

Algebraic multigrid (AMG) solvers and preconditioners are some of the fastest numerical methods to solve linear systems, particularly in a parallel environment, scaling to hundreds of thousands of cores. Most AMG methods and theory assume a symmetric positive definite operator. This paper presents a new variation on classical AMG for nonsymmetric matrices (denoted lAIR), based on a local approximation to the ideal restriction operator, coupled with F-relaxation. A new block decomposition of the AMG error-propagation operator is used for a spectral analysis of convergence, and the efficacy of the algorithm is demonstrated on systems arising from the discrete form of the advection-diffusion-reaction equation. lAIR is shown to be a robust solver for various discretizations of the advection-diffusion-reaction equation, including time-dependent and steady-state, from purely advective to purely diffusive. Convergence is robust for discretizations on unstructured meshes and using higher-order finite elements, and is particularly effective on upwind discontinuous Galerkin discretizations. Although the implementation used here is not parallel, each part of the algorithm is highly parallelizable, avoiding common multigrid adjustments for strong advection such as line-relaxation and K- or W-cycles that can be effective in serial, but suffer from high communication costs in parallel, limiting their scalability.

Quasi-Monte Carlo (QMC) method is a useful numerical tool for pricing and hedging of complex financial derivatives. These problems are usually of high dimensionality and discontinuities. The two factors may significantly deteriorate the performance of the QMC method. This paper develops an integrated method that overcomes the challenges of the high dimensionality and discontinuities concurrently. For this purpose, a smoothing method is proposed to remove the discontinuities for some typical functions arising from financial engineering. To make the smoothing method applicable for more general functions, a new path generation method is designed for simulating the paths of the underlying assets such that the resulting function has the required form. The new path generation method has an additional power to reduce the effective dimension of the target function. Our proposed method caters for a large variety of model specifications, including the Black-Scholes, exponential normal inverse Gaussian L\'evy, and Heston models. Numerical experiments dealing with these models show that in the QMC setting the proposed smoothing method in combination with the new path generation method can lead to a dramatic variance reduction for pricing exotic options with discontinuous payoffs and for calculating options' Greeks. The investigation on the effective dimension and the related characteristics explains the significant enhancement of the combined procedure.

One way of getting insight into non-Gaussian measures, posed on infinite dimensional Hilbert spaces, is to first obtain best fit Gaussian approximations, which are more amenable to numerical approximation. These Gaussians can then be used to accelerate sampling algorithms. This begs the questions of how one should measure optimality and how the optimizers can be obtained. Here, we consider the problem of minimizing the distance with respect to relative entropy. We examine this minimization problem by seeking roots of the first variation of relative entropy, taken with respect to the mean of the Gaussian, leaving the covariance fixed. Adapting a convergence analysis of Robbins-Monro to the infinite dimensional setting, we can justify the application of this algorithm and highlight necessary assumptions to ensure convergence, not only in the context of relative entropy minimization, but other infinite dimensional problems as well. Numerical examples in path space, showing the robustness of this method with respect to dimension, are provided.

The conditioning of implicit Runge-Kutta (RK) integration for linear finite element approximation of diffusion equations on general anisotropic meshes is investigated. Bounds are established for the condition number of the resulting linear system with and without diagonal preconditioning for the implicit Euler and general implicit RK methods. Two solution strategies are considered for the linear system resulting from general implicit RK integration: the simultaneous solution (the system is solved as a whole) and a successive solution which follows the commonly used implementation of implicit RK methods to first transform the system into smaller systems using the Jordan normal form of the RK matrix and then solve them successively. For the simultaneous solution in case of a positive semidefinite symmetric part of the RK coefficient matrix and for the successive solution it is shown that . If the smallest eigenvalue of the symmetric part of the RK coefficient matrix is negative and the simultaneous solution strategy is used, an upper bound on the time step is given so that the system matrix is positive definite. The obtained bounds for the condition number have explicit geometric interpretations and take the interplay between the diffusion matrix and the mesh geometry into full consideration. They show that there are three mesh-dependent factors that can affect the conditioning: the number of elements, the mesh nonuniformity measured in the Euclidean metric, and the mesh nonuniformity with respect to the inverse of the diffusion matrix. They also reveal that the preconditioning using the diagonal of the system matrix, the mass matrix, or the lumped mass matrix can effectively eliminate the effects of the mesh nonuniformity measured in the Euclidean metric. Numerical examples are given.

We consider solving the exterior Dirichlet problem for the Helmholtz equation with the $h$-version of the boundary element method (BEM) using the standard second-kind combined-field integral equations. We prove a new, sharp bound on how the number of GMRES iterations must grow with the wavenumber $k$ to have the error in the iterative solution bounded independently of $k$ as $k\rightarrow \infty$ when the boundary of the obstacle is analytic and has strictly positive curvature. To our knowledge, this result is the first-ever sharp bound on how the number of GMRES iterations depends on the wavenumber for an integral equation used to solve a scattering problem. We also prove new bounds on how $h$ must decrease with $k$ to maintain $k$-independent quasi-optimality of the Galerkin solutions as $k \rightarrow \infty$ when the obstacle is nontrapping.

We introduce a new family of numerical algorithms for approximating solutions of general high-dimensional semilinear parabolic partial differential equations at single space-time points. The algorithm is obtained through a delicate combination of the Feynman-Kac and the Bismut-Elworthy-Li formulas, and an approximate decomposition of the Picard fixed-point iteration with multilevel accuracy. The algorithm has been tested on a variety of semilinear partial differential equations that arise in physics and finance, with very satisfactory results. Analytical tools needed for the analysis of such algorithms, including a semilinear Feynman-Kac formula, a new class of semi-norms and their recursive inequalities, are also introduced. They allow us to prove for semilinear heat equations with gradient-independent nonlinearity that the computational complexity of the proposed algorithm is bounded by $O(d\,\varepsilon^{-(4+\delta)})$ for any $\delta \in (0,\infty)$ under suitable assumptions, where $d\in \mathbb{N}$ is the dimensionality of the problem and $\varepsilon\in(0,\infty)$ is the prescribed accuracy.

Temporal or spatial structures are readily extracted from complex data by modal decompositions like Proper Orthogonal Decomposition (POD) or Dynamic Mode Decomposition (DMD). Subspaces of such decompositions serve as reduced order models and define either spatial structures in time or temporal structures in space. On the contrary, convecting phenomena pose a major problem to those decompositions. A structure traveling with a certain group velocity will be perceived as a plethora of modes in time or space respectively. This manifests itself for example in poorly decaying singular values when using a POD. The poor decay is counter-intuitive, since a single structure is expected to be represented by a few modes. The intuition proves to be correct and we show that in a properly chosen reference frame along the characteristics defined by the group velocity, a POD or DMD reduces moving structures to a few modes, as expected. Beyond serving as a reduced model, the resulting entity can be used to define a constant or minimally changing structure in turbulent flows. This can be interpreted as an empirical counterpart to exact coherent structures. We present the method and its application to a head vortex of a compressible starting jet.

In this work, we study the numerical approximation of local fluctuations of certain classes of parabolic stochastic partial differential equations (SPDEs). Our focus is on effects for small spatially-correlated noise on a time scale before large deviation effects have occurred. In particular, we are interested in the local directions of the noise described by a covariance operator. We introduce a new strategy and prove a Combined ERror EStimate (CERES) for the five main errors: the spatial discretization error, the local linearization error, the noise truncation error, the local relaxation error to steady state, and the approximation error via an iterative low-rank matrix algorithm. In summary, we obtain one CERES describing, apart from modelling of the original equations and standard round-off, all sources of error for a local fluctuation analysis of an SPDE in one estimate. To prove our results, we rely on a combination of methods from optimal Galerkin approximation of SPDEs, covariance moment estimates, analytical techniques for Lyapunov equations, iterative numerical schemes for low-rank solution of Lyapunov equations, and working with related spectral norms for different classes of operators.

In this paper we present algorithms for an efficient implementation of the Localized Orthogonal Decomposition method (LOD). The LOD is a multiscale method for the numerical simulation of partial differential equations with a continuum of inseparable scales. We show how the method can be implemented in a fairly standard Finite Element framework and discuss its realization for different types of problems, such as linear elliptic problems with rough coefficients and linear eigenvalue problems.

We develop a framework that allows the use of the multi-level Monte Carlo (MLMC) methodology (Giles2015) to calculate expectations with respect to the invariant measure of an ergodic SDE. In that context, we study the (over-damped) Langevin equations with a strongly concave potential. We show that, when appropriate contracting couplings for the numerical integrators are available, one can obtain a uniform in time estimate of the MLMC variance in contrast to the majority of the results in the MLMC literature. As a consequence, a root mean square error of $\mathcal{O}(\varepsilon)$ is achieved with $\mathcal{O}(\varepsilon^{-2})$ complexity on par with Markov Chain Monte Carlo (MCMC) methods, which however can be computationally intensive when applied to large data sets. Finally, we present a multi-level version of the recently introduced Stochastic Gradient Langevin Dynamics (SGLD) method (Welling and Teh, 2011) built for large datasets applications. We show that this is the first stochastic gradient MCMC method with complexity $\mathcal{O}(\varepsilon^{-2}|\log {\varepsilon}|^{3})$, in contrast to the complexity $\mathcal{O}(\varepsilon^{-3})$ of currently available methods. Numerical experiments confirm our theoretical findings.

Deep neural networks have become invaluable tools for supervised machine learning, e.g., classification of text or images. While often offering superior results over traditional techniques and successfully expressing complicated patterns in data, deep architectures are known to be challenging to design and train such that they generalize well to new data. Important issues with deep architectures are numerical instabilities in derivative-based learning algorithms commonly called exploding or vanishing gradients. In this paper we propose new forward propagation techniques inspired by systems of Ordinary Differential Equations (ODE) that overcome this challenge and lead to well-posed learning problems for arbitrarily deep networks. The backbone of our approach is our interpretation of deep learning as a parameter estimation problem of nonlinear dynamical systems. Given this formulation, we analyze stability and well-posedness of deep learning and use this new understanding to develop new network architectures. We relate the exploding and vanishing gradient phenomenon to the stability of the discrete ODE and present several strategies for stabilizing deep learning for very deep networks. While our new architectures restrict the solution space, several numerical experiments show their competitiveness with state-of-the-art networks.

The generalized Debye source representation of time-harmonic electromagnetic fields yields well-conditioned second-kind integral equations for a variety of boundary value problems, including the problems of scattering from perfect electric conductors and dielectric bodies. Furthermore, these representations, and resulting integral equations, are fully stable in the static limit as $\omega \to 0$ in multiply connected geometries. In this paper, we present the first high-order accurate solver based on this representation for bodies of revolution. The resulting solver uses a Nystr\"om discretization of a one-dimensional generating curve and high-order integral equation methods for applying and inverting surface differentials. The accuracy and speed of the solvers are demonstrated in several numerical examples.

In this paper we present criteria for the choice of the shape parameter c contained in the famous radial function multiquadric. It may be of interest to RBF people and all people using radial basis functions to do approximation.

The topic of these notes could be easily expanded into a full one-semester course. Nevertheless, we shall try to give some flavour along with theoretical bases of spectral and pseudo-spectral methods. The main focus is made on Fourier-type discretizations, even if some indications on how to handle non-periodic problems via Tchebyshev and Legendre approaches are made as well. The applications presented here are diffusion-type problems in accordance with the topics of the PhD school.

In this paper, we prove necessary and sufficient conditions for a hybridizable discontinuous Galerkin (HDG) method to satisfy a multisymplectic conservation law, when applied to a canonical Hamiltonian system of partial differential equations. We show that these conditions are satisfied by the "hybridized" versions of several of the most commonly-used finite element methods, including mixed, nonconforming, and discontinuous Galerkin methods. (Interestingly, for the continuous Galerkin method in dimension greater than one, we show that multisymplecticity only holds in a weaker sense.) Consequently, these general-purpose finite element methods may be used for structure-preserving discretization (or semidiscretization) of canonical Hamiltonian systems of ODEs or PDEs. This establishes multisymplecticity for a large class of arbitrarily-high-order methods on unstructured meshes.