scispace - formally typeset
Search or ask a question

Showing papers in "arXiv: Numerical Analysis in 2019"


Posted Content
TL;DR: Finite element methods for approximating the time harmonic Maxwell equations and Discontinuous Galerkin methods are surveyed, focusing on comparing error estimates for problems with spatially varying coefficients.
Abstract: We survey finite element methods for approximating the time harmonic Maxwell equations. We concentrate on comparing error estimates for problems with spatially varying coefficients. For the conforming edge finite element methods, such estimates allow, at least, piecewise smooth coefficients. But for Discontinuous Galerkin (DG) methods, the state of the art of error analysis is less advanced (we consider three DG families of methods: Interior Penalty type, Hybridizable DG, and Trefftz type methods). Nevertheless, DG methods offer significant potential advantages compared to conforming methods.

473 citations


Journal ArticleDOI
TL;DR: The approximation power of deep feed-forward neural networks is quantitatively characterizes in terms of the number of neurons, i.e., the product of the network width and depth, to provide a general guide for selecting the width and the depth of ReLU FNNs to approximate continuous functions especially in parallel computing.
Abstract: This paper quantitatively characterizes the approximation power of deep feed-forward neural networks (FNNs) in terms of the number of neurons. It is shown by construction that ReLU FNNs with width $\mathcal{O}\big(\max\{d\lfloor N^{1/d}\rfloor,\, N+1\}\big)$ and depth $\mathcal{O}(L)$ can approximate an arbitrary Holder continuous function of order $\alpha\in (0,1]$ on $[0,1]^d$ with a nearly tight approximation rate $\mathcal{O}\big(\sqrt{d} N^{-2\alpha/d}L^{-2\alpha/d}\big)$ measured in $L^p$-norm for any $N,L\in \mathbb{N}^+$ and $p\in[1,\infty]$. More generally for an arbitrary continuous function $f$ on $[0,1]^d$ with a modulus of continuity $\omega_f(\cdot)$, the constructive approximation rate is $\mathcal{O}\big(\sqrt{d}\,\omega_f( N^{-2/d}L^{-2/d})\big)$. We also extend our analysis to $f$ on irregular domains or those localized in an $\varepsilon$-neighborhood of a $d_{\mathcal{M}}$-dimensional smooth manifold $\mathcal{M}\subseteq [0,1]^d$ with $d_{\mathcal{M}}\ll d$. Especially, in the case of an essentially low-dimensional domain, we show an approximation rate $\mathcal{O}\big(\omega_f(\tfrac{\varepsilon}{1-\delta}\sqrt{\tfrac{d}{d_\delta}}+\varepsilon)+\sqrt{d}\,\omega_f(\tfrac{\sqrt{d}}{(1-\delta)\sqrt{d_\delta}}N^{-2/d_\delta}L^{-2/d_\delta})\big)$ for ReLU FNNs to approximate $f$ in the $\varepsilon$-neighborhood, where $d_\delta=\mathcal{O}\big(d_{\mathcal{M}}\tfrac{\ln (d/\delta)}{\delta^2}\big)$ for any $\delta\in(0,1)$ as a relative error for a projection to approximate an isometry when projecting $\mathcal{M}$ to a $d_{\delta}$-dimensional domain.

101 citations


Journal ArticleDOI
TL;DR: A numerical method that divides the PDE approximation problem into a sequence of separate learning problems that combines operator splitting with deep learning and can handle extremely high-dimensional PDEs.
Abstract: In this paper we introduce a numerical method for nonlinear parabolic PDEs that combines operator splitting with deep learning. It divides the PDE approximation problem into a sequence of separate learning problems. Since the computational graph for each of the subproblems is comparatively small, the approach can handle extremely high-dimensional PDEs. We test the method on different examples from physics, stochastic control and mathematical finance. In all cases, it yields very good results in up to 10,000 dimensions with short run times.

94 citations


Posted Content
TL;DR: In this article, the authors derive upper bounds on the complexity of ReLU neural networks approximating the solution maps of parametric partial differential equations, without any knowledge of its concrete shape, using the inherent low-dimensionality of the solution manifold to obtain approximation rates which are significantly superior to those provided by classical neural networks.
Abstract: We derive upper bounds on the complexity of ReLU neural networks approximating the solution maps of parametric partial differential equations. In particular, without any knowledge of its concrete shape, we use the inherent low-dimensionality of the solution manifold to obtain approximation rates which are significantly superior to those provided by classical neural network approximation results. Concretely, we use the existence of a small reduced basis to construct, for a large variety of parametric partial differential equations, neural networks that yield approximations of the parametric solution maps in such a way that the sizes of these networks essentially only depend on the size of the reduced basis.

77 citations


Journal ArticleDOI
TL;DR: The numerical results clearly show that the Deep Nitsche Method is naturally nonlinear, naturally adaptive and has the potential to work on rather high dimensions.
Abstract: We propose a new method to deal with the essential boundary conditions encountered in the deep learning-based numerical solvers for partial differential equations. The trial functions representing by deep neural networks are non-interpolatory, which makes the enforcement of the essential boundary conditions a nontrivial matter. Our method resorts to Nitsche's variational formulation to deal with this difficulty, which is consistent, and does not require significant extra computational costs. We prove the error estimate in the energy norm and illustrate the method on several representative problems posed in at most 100 dimension.

72 citations


Posted Content
TL;DR: A fast recursive algorithm is provided for efficient evaluation of the derivation of polynomial basis functions for approximating isometry and permutation invariant functions, particularly with an eye to modelling properties of atomistic systems.
Abstract: The Atomic Cluster Expansion (Drautz, Phys. Rev. B 99, 2019) provides a framework to systematically derive polynomial basis functions for approximating isometry and permutation invariant functions, particularly with an eye to modelling properties of atomistic systems. Our presentation extends the derivation by proposing a precomputation algorithm that yields immediate guarantees that a complete basis is obtained. We provide a fast recursive algorithm for efficient evaluation and illustrate its performance in numerical tests. Finally, we discuss generalisations and open challenges, particularly from a numerical stability perspective, around basis optimisation and parameter estimation, paving the way towards a comprehensive analysis of the convergence to a high-fidelity reference model.

67 citations


Posted Content
TL;DR: In this paper, a third order accurate exponential time differencing (ETD) numerical scheme for the no-slope-selection (NSS) equation of the epitaxial thin film growth model, with Fourier pseudo-spectral discretization in space, was proposed and analyzed.
Abstract: In this paper we propose and analyze a (temporally) third order accurate exponential time differencing (ETD) numerical scheme for the no-slope-selection (NSS) equation of the epitaxial thin film growth model, with Fourier pseudo-spectral discretization in space. A linear splitting is applied to the physical model, and an ETD-based multistep approximation is used for time integration of the corresponding equation. In addition, a third order accurate Douglas-Dupont regularization term, in the form of $-A \dt^2 \phi_0 (L_N) \Delta_N^2 ( u^{n+1} - u^n)$, is added in the numerical scheme. A careful Fourier eigenvalue analysis results in the energy stability in a modified version, and a theoretical justification of the coefficient $A$ becomes available. As a result of this energy stability analysis, a uniform in time bound of the numerical energy is obtained. And also, the optimal rate convergence analysis and error estimate are derived in details, in the $\ell^\infty (0,T; H_h^1) \cap \ell^2 (0,T; H_h^3)$ norm, with the help of a careful eigenvalue bound estimate, combined with the nonlinear analysis for the NSS model. This convergence estimate is the first such result for a third order accurate scheme for a gradient flow. Some numerical simulation results are presented to demonstrate the efficiency of the numerical scheme and the third order convergence. The long time simulation results for $\varepsilon=0.02$ (up to $T=3 \times 10^5$) have indicated a logarithm law for the energy decay, as well as the power laws for growth of the surface roughness and the mound width. In particular, the power index for the surface roughness and the mound width growth, created by the third order numerical scheme, is more accurate than those produced by certain second order energy stable schemes in the existing literature.

66 citations


Journal ArticleDOI
TL;DR: A continuous formulation of machine learning, as a problem in the calculus of variations and differential-integral equations, in the spirit of classical numerical analysis is presented.
Abstract: We present a continuous formulation of machine learning, as a problem in the calculus of variations and differential-integral equations, in the spirit of classical numerical analysis. We demonstrate that conventional machine learning models and algorithms, such as the random feature model, the two-layer neural network model and the residual neural network model, can all be recovered (in a scaled form) as particular discretizations of different continuous formulations. We also present examples of new models, such as the flow-based random feature model, and new algorithms, such as the smoothed particle method and spectral method, that arise naturally from this continuous formulation. We discuss how the issues of generalization error and implicit regularization can be studied under this framework.

55 citations


Posted Content
TL;DR: A theorem concerning the approximation of bandlimited multivariate functions by deep ReLU networks for which the curse of the dimensionality is overcome is proved.
Abstract: We prove a theorem concerning the approximation of bandlimited multivariate functions by deep ReLU networks for which the curse of the dimensionality is overcome. Our theorem is based on a result by Maurey and on the ability of deep ReLU networks to approximate Chebyshev polynomials and analytic functions efficiently.

51 citations


Posted Content
TL;DR: Yan et al. as discussed by the authors developed a new numerical method for nonlinear, nonlocal partial differential equations, arising in models of porous media, materials science, and biological swarming, combining the classical theory of optimal transport with modern operator splitting techniques.
Abstract: Combining the classical theory of optimal transport with modern operator splitting techniques, we develop a new numerical method for nonlinear, nonlocal partial differential equations, arising in models of porous media, materials science, and biological swarming. Our method proceeds as follows: First, we discretize in time, either via the classical JKO scheme or via a novel Crank-Nicolson type method we introduce. Next, we use the Benamou-Brenier dynamical characterization of the Wasserstein distance to reduce computing the solution of the discrete time equations to solving fully discrete minimization problems, with strictly convex objective functions and linear constraints. Third, we compute the minimizers by applying a recently introduced, provably convergent primal dual splitting scheme for three operators [Yan 2018]. By leveraging the PDEs' underlying variational structure, our method overcomes stability issues present in previous numerical work built on explicit time discretizations, which suffer due to the equations' strong nonlinearities and degeneracies. Our method is also naturally positivity and mass preserving and, in the case of the JKO scheme, energy decreasing. We prove that minimizers of the fully discrete problem converge to minimizers of the spatially continuous, discrete time problem as the spatial discretization is refined. We conclude with simulations of nonlinear PDEs and Wasserstein geodesics in one and two dimensions that illustrate the key properties of our approach, including higher order convergence our novel Crank-Nicolson type method, when compared to the classical JKO method.

49 citations


Posted Content
TL;DR: In this article, it was shown that for a wide class of controlled stochastic differential equations (SDEs) with stiff coefficients, the value functions of corresponding zero-sum games can be represented by a deep artificial neural network (DNN), whose complexity grows at most polynomially in both the dimension of the state equation and the reciprocal of the required accuracy.
Abstract: In this paper, we establish that for a wide class of controlled stochastic differential equations (SDEs) with stiff coefficients, the value functions of corresponding zero-sum games can be represented by a deep artificial neural network (DNN), whose complexity grows at most polynomially in both the dimension of the state equation and the reciprocal of the required accuracy. Such nonlinear stiff systems may arise, for example, from Galerkin approximations of controlled stochastic partial differential equations (SPDEs), or controlled PDEs with uncertain initial conditions and source terms. This implies that DNNs can break the curse of dimensionality in numerical approximations and optimal control of PDEs and SPDEs. The main ingredient of our proof is to construct a suitable discrete-time system to effectively approximate the evolution of the underlying stochastic dynamics. Similar ideas can also be applied to obtain expression rates of DNNs for value functions induced by stiff systems with regime switching coefficients and driven by general Levy noise.

Posted Content
TL;DR: In this article, the authors propose a framework for learning multigrid solvers for solving large-scale partial differential equations (PDEs) by learning a mapping from a family of parameterized PDEs to prolongation operators using an efficient and unsupervised loss function.
Abstract: Constructing fast numerical solvers for partial differential equations (PDEs) is crucial for many scientific disciplines. A leading technique for solving large-scale PDEs is using multigrid methods. At the core of a multigrid solver is the prolongation matrix, which relates between different scales of the problem. This matrix is strongly problem-dependent, and its optimal construction is critical to the efficiency of the solver. In practice, however, devising multigrid algorithms for new problems often poses formidable challenges. In this paper we propose a framework for learning multigrid solvers. Our method learns a (single) mapping from a family of parameterized PDEs to prolongation operators. We train a neural network once for the entire class of PDEs, using an efficient and unsupervised loss function. Experiments on a broad class of 2D diffusion problems demonstrate improved convergence rates compared to the widely used Black-Box multigrid scheme, suggesting that our method successfully learned rules for constructing prolongation matrices.

Posted Content
TL;DR: The subject of the main result of this article is to provide space-time error estimates for DNN approximation of Euler approximations of certain perturbed differential equations.
Abstract: Over the last few years deep artificial neural networks (DNNs) have very successfully been used in numerical simulations for a wide variety of computational problems including computer vision, image classification, speech recognition, natural language processing, as well as computational advertisement. In addition, it has recently been proposed to approximate solutions of partial differential equations (PDEs) by means of stochastic learning problems involving DNNs. There are now also a few rigorous mathematical results in the scientific literature which provide error estimates for such deep learning based approximation methods for PDEs. All of these articles provide spatial error estimates for neural network approximations for PDEs but do not provide error estimates for the entire space-time error for the considered neural network approximations. It is the subject of the main result of this article to provide space-time error estimates for DNN approximations of Euler approximations of certain perturbed differential equations. Our proof of this result is based (i) on a certain artificial neural network (ANN) calculus and (ii) on ANN approximation results for products of the form $[0,T]\times \mathbb{R}^d i (t,x)\mapsto tx\in \mathbb{R}^d$ where $T\in (0,\infty)$, $d\in \mathbb{N}$, which we both develop within this article.

Posted Content
TL;DR: The Kolmogorov N -width d N ( M ) describes the rate of the worst-case error arising from a projection onto the best-possible linear subspace of H of dimension N ∈ N and sets a limit to any projection-based approximation such as determined by the reduced basis method.
Abstract: The Kolmogorov $N$-width $d_N(\mathcal{M})$ describes the rate of the worst-case error (w.r.t.\ a subset $\mathcal{M}\subset H$ of a normed space $H$) arising from a projection onto the best-possible linear subspace of $H$ of dimension $N\in\mathbb{N}$. Thus, $d_N(\mathcal{M})$ sets a limit to any projection-based approximation such as determined by the reduced basis method. While it is known that $d_N(\mathcal{M})$ decays exponentially fast for many linear coercive parametrized partial differential equations, i.e., $d_N(\mathcal{M})=\mathcal{O}(e^{-\beta N})$, we show in this note, that only $d_N(\mathcal{M}) =\mathcal{O}(N^{-1/2})$ for initial-boundary-value problems of the hyperbolic wave equation with discontinuous initial conditions. This is aligned with the known slow decay of $d_N(\mathcal{M})$ for the linear transport problem.

Journal ArticleDOI
TL;DR: An energy stable numerical scheme for the square phase field crystal (SPFC) equation, a gradient flow modeling crystal dynamics at the atomic scale in space but on diffusive scales in time, with unique solvability, energy stability, and an optimal rate convergence analysis is derived.
Abstract: In this paper we propose and analyze an energy stable numerical scheme for the square phase field crystal (SPFC) equation, a gradient flow modeling crystal dynamics at the atomic scale in space but on diffusive scales in time. In particular, a modification of the free energy potential to the standard phase field crystal model leads to a composition of the 4-Laplacian and the regular Laplacian operators. To overcome the difficulties associated with this highly nonlinear operator, we design numerical algorithms based on the structures of the individual energy terms. A Fourier pseudo-spectral approximation is taken in space, in such a way that the energy structure is respected, and summation-by-parts formulae enable us to study the discrete energy stability for such a high-order spatial discretization. In the temporal approximation, a second order BDF stencil is applied, combined with an appropriate extrapolation for the concave diffusion term(s). A second order artificial Douglas-Dupont-type regularization term is added to ensure energy stability, and a careful analysis leads to the artificial linear diffusion coming at an order lower that that of surface diffusion term. Such a choice leads to reduced numerical dissipation. At a theoretical level, the unique solvability, energy stability are established, and an optimal rate convergence analysis is derived in the $\ell^\infty (0,T; \ell^2) \cap \ell^2 (0,T; H_N^3)$ norm. In the numerical implementation, the preconditioned steepest descent (PSD) iteration is applied to solve for the composition of the highly nonlinear 4-Laplacian term and the standard Laplacian term, and a geometric convergence is assured for such an iteration. Finally, a few numerical experiments are presented, which confirm the robustness and accuracy of the proposed scheme.

Posted Content
TL;DR: This work demonstrates how further regularization can be imposed, incorporating prior information about the underlying unknown, and studies how to impose Tikhonov-like Sobolev penalties in the ensemble inversion context.
Abstract: Ensemble Kalman inversion is a parallelizable methodology for solving inverse or parameter estimation problems. Although it is based on ideas from Kalman filtering, it may be viewed as a derivative-free optimization method. In its most basic form it regularizes ill-posed inverse problems through the subspace property: the solution found is in the linear span of the initial ensemble employed. In this work we demonstrate how further regularization can be imposed, incorporating prior information about the underlying unknown. In particular we study how to impose Tikhonov-like Sobolev penalties. As well as introducing this modified ensemble Kalman inversion methodology, we also study its continuous-time limit, proving ensemble collapse; in the language of multi-agent optimization this may be viewed as reaching consensus. We also conduct a suite of numerical experiments to highlight the benefits of Tikhonov regularization in the ensemble inversion context.

Posted Content
TL;DR: It is found that the block circulant operator established an isomorphism between tensors and matrices that is used to prove the F-stochastic structure is invariant under generalized tensor functions.
Abstract: In this paper, we present the definition of generalized tensor function according to the tensor singular value decomposition (T-SVD) via the tensor T-product. Also, we introduce the compact singular value decomposition (T-CSVD) of tensors via the T-product, from which the projection operators and Moore Penrose inverse of tensors are also obtained. We also establish the Cauchy integral formula for tensors by using the partial isometry tensors and applied it into the solution of tensor equations. Then we establish the generalized tensor power and the Taylor expansion of tensors. Explicit generalized tensor functions are also listed. We define the tensor bilinear and sesquilinear forms and proposed theorems on structures preserved by generalized tensor functions. For complex tensors, we established an isomorphism between complex tensors and real tensors. In the last part of our paper, we find that the block circulant operator established an isomorphism between tensors and matrices. This isomorphism is used to prove the F-stochastic structure is invariant under generalized tensor functions. The concept of invariant tensor cones is also raised.

Posted Content
TL;DR: From the recovery analysis of Lasso, a new definition of Noise-to-Signal ratio is proposed, which better represents the level of noise in the case of PDE identification, and a performance guarantee is established based on an incoherence property.
Abstract: Identifying unknown differential equations from a given set of discrete time dependent data is a challenging problem. A small amount of noise can make the recovery unstable, and nonlinearity and differential equations with varying coefficients add complexity to the problem. We assume that the governing partial differential equation (PDE) is a linear combination of a subset of a prescribed dictionary containing different differential terms, and the objective of this paper is to find the correct coefficients. We propose a new direction based on the fundamental idea of convergence analysis of numerical PDE schemes. We utilize Lasso for efficiency, and a performance guarantee is established based on an incoherence property. The main contribution is to validate and correct the results by Time Evolution Error (TEE). The new algorithm, called Identifying Differential Equations with Numerical Time evolution (IDENT), is explored for data with non-periodic boundary conditions, noisy data and PDEs with varying coefficients. From the recovery analysis of Lasso, we propose a new definition of Noise-to-Signal ratio, which better represents the level of noise in the case of PDE identification. We systematically analyze the effects of data generations and downsampling, and propose an order preserving denoising method called Least-Squares Moving Average (LSMA), to preprocess the given data. For the identification of PDEs with varying coefficients, we propose to add Base Element Expansion (BEE) to aide the computation. Various numerical experiments from basic tests to noisy data, downsampling effects and varying coefficients are presented.

Posted Content
TL;DR: Finite element results obtained reveal the robustness of quasi-Newton monolithic schemes, with convergence readily attained under both stable and unstable cracking conditions, and since the solution method is unconditionally stable, very significant computational gains are observed relative to the widely used staggered solution schemes.
Abstract: We investigate the potential of quasi-Newton methods in facilitating convergence of monolithic solution schemes for phase field fracture modelling. Several paradigmatic boundary value problems are addressed, spanning the fields of quasi-static fracture, fatigue damage and dynamic cracking. The finite element results obtained reveal the robustness of quasi-Newton monolithic schemes, with convergence readily attained under both stable and unstable cracking conditions. Moreover, since the solution method is unconditionally stable, very significant computational gains are observed relative to the widely used staggered solution schemes. In addition, a new adaptive time increment scheme is presented to further reduces the computational cost while allowing to accurately resolve sudden changes in material behavior, such as unstable crack growth. Computation times can be reduced by several orders of magnitude, with the number of load increments required by the corresponding staggered solution being up to 3000 times higher. Quasi-Newton monolithic solution schemes can be a key enabler for large scale phase field fracture simulations. Implications are particularly relevant for the emerging field of phase field fatigue, as results show that staggered cycle-by-cycle calculations are prohibitive in mid or high cycle fatigue. The finite element codes are available to download from this http URL.

Journal ArticleDOI
TL;DR: It is proved for the first time that in the case of semilinear heat equations with gradient-independent nonlinearities that the numbers of parameters of the employed deep neural networks grow at most polynomially in both the PDE dimension and the reciprocal of the prescribed approximation accuracy.
Abstract: Deep neural networks and other deep learning methods have very successfully been applied to the numerical approximation of high-dimensional nonlinear parabolic partial differential equations (PDEs), which are widely used in finance, engineering, and natural sciences. In particular, simulations indicate that algorithms based on deep learning overcome the curse of dimensionality in the numerical approximation of solutions of semilinear PDEs. For certain linear PDEs this has also been proved mathematically. The key contribution of this article is to rigorously prove this for the first time for a class of nonlinear PDEs. More precisely, we prove in the case of semilinear heat equations with gradient-independent nonlinearities that the numbers of parameters of the employed deep neural networks grow at most polynomially in both the PDE dimension and the reciprocal of the prescribed approximation accuracy. Our proof relies on recently introduced multilevel Picard approximations of semilinear PDEs.

Journal ArticleDOI
TL;DR: It is shown that deep networks with rectified power units (RePU) can give better approximations for smooth functions than deep ReLU networks, and some efficient algorithms proposed in this paper to convert polynomials into deep RePU networks of optimal size with no approximation error are proposed.
Abstract: Deep neural networks with rectified linear units (ReLU) are getting more and more popular due to their universal representation power and successful applications. Some theoretical progress regarding the approximation power of deep ReLU network for functions in Sobolev space and Korobov space have recently been made by [D. Yarotsky, Neural Network, 94:103-114, 2017] and [H. Montanelli and Q. Du, SIAM J Math. Data Sci., 1:78-92, 2019], etc. In this paper, we show that deep networks with rectified power units (RePU) can give better approximations for smooth functions than deep ReLU networks. Our analysis bases on classical polynomial approximation theory and some efficient algorithms proposed in this paper to convert polynomials into deep RePU networks of optimal size with no approximation error. Comparing to the results on ReLU networks, the sizes of RePU networks required to approximate functions in Sobolev space and Korobov space with an error tolerance $\varepsilon$, by our constructive proofs, are in general $\mathcal{O}(\log\frac{1}{\varepsilon})$ times smaller than the sizes of corresponding ReLU networks constructed in most of the existing literature. Comparing to the classical results of Mhaskar [Mhaskar, Adv. Comput. Math. 1:61-80, 1993], our constructions use less number of activation functions and numerically more stable, they can be served as good initials of deep RePU networks and further trained to break the limit of linear approximation theory. The functions represented by RePU networks are smooth functions, so they naturally fit in the places where derivatives are involved in the loss function.

Journal ArticleDOI
TL;DR: In this paper, it was shown that the problem of solving semilinear Black-Scholes equations is a polynomially tractable approximation problem and that the computational effort of the algorithm grows polynomial both in the dimension and the reciprocal of the prescribed approximation accuracy.
Abstract: Parabolic partial differential equations (PDEs) are widely used in the mathematical modeling of natural phenomena and man made complex systems. In particular, parabolic PDEs are a fundamental tool to determine fair prices of financial derivatives in the financial industry. The PDEs appearing in financial engineering applications are often nonlinear and high dimensional since the dimension typically corresponds to the number of considered financial assets. A major issue is that most approximation methods for nonlinear PDEs in the literature suffer under the so-called curse of dimensionality in the sense that the computational effort to compute an approximation with a prescribed accuracy grows exponentially in the dimension of the PDE or in the reciprocal of the prescribed approximation accuracy and nearly all approximation methods have not been shown not to suffer under the curse of dimensionality. Recently, a new class of approximation schemes for semilinear parabolic PDEs, termed full history recursive multilevel Picard (MLP) algorithms, were introduced and it was proven that MLP algorithms do overcome the curse of dimensionality for semilinear heat equations. In this paper we extend those findings to a more general class of semilinear PDEs including as special cases semilinear Black-Scholes equations used for the pricing of financial derivatives with default risks. More specifically, we introduce an MLP algorithm for the approximation of solutions of semilinear Black-Scholes equations and prove that the computational effort of our method grows at most polynomially both in the dimension and the reciprocal of the prescribed approximation accuracy. This is, to the best of our knowledge, the first result showing that the approximation of solutions of semilinear Black-Scholes equations is a polynomially tractable approximation problem.

Journal ArticleDOI
TL;DR: In this article, a phase-field method for continuous modeling of cracks with frictional contacts is proposed, which can represent arbitrary crack geometry without an explicit function or basis enrichment, and does not require an algorithm for imposing contact constraints.
Abstract: We introduce a phase-field method for continuous modeling of cracks with frictional contacts. Compared with standard discrete methods for frictional contacts, the phase-field method has two attractive features: (1) it can represent arbitrary crack geometry without an explicit function or basis enrichment, and (2) it does not require an algorithm for imposing contact constraints. The first feature, which is common in phase-field models of fracture, is attained by regularizing a sharp interface geometry using a surface density functional. The second feature, which is a unique advantage for contact problems, is achieved by a new approach that calculates the stress tensor in the regularized interface region depending on the contact condition of the interface. Particularly, under a slip condition, this approach updates stress components in the slip direction using a standard contact constitutive law, while making other stress components compatible with stress in the bulk region to ensure non-penetrating deformation in other directions. We verify the proposed phase-field method using stationary interface problems simulated by discrete methods in the literature. Subsequently, by allowing the phase field to evolve according to brittle fracture theory, we demonstrate the proposed method's capability for modeling crack growth with frictional contact.

Posted Content
TL;DR: This work presents randomized versions of two well-known compression algorithms, namely, HOSVD and STHOSVD, and develops variants of these algorithms that tackle specific challenges posed by large-scale datasets.
Abstract: Many applications in data science and scientific computing involve large-scale datasets that are expensive to store and compute with, but can be efficiently compressed and stored in an appropriate tensor format. In recent years, randomized matrix methods have been used to efficiently and accurately compute low-rank matrix decompositions. Motivated by this success, we focus on developing randomized algorithms for tensor decompositions in the Tucker representation. Specifically, we present randomized versions of two well-known compression algorithms, namely, HOSVD and STHOSVD. We present a detailed probabilistic analysis of the error of the randomized tensor algorithms. We also develop variants of these algorithms that tackle specific challenges posed by large-scale datasets. The first variant adaptively finds a low-rank representation satisfying a given tolerance and it is beneficial when the target-rank is not known in advance. The second variant preserves the structure of the original tensor, and is beneficial for large sparse tensors that are difficult to load in memory. We consider several different datasets for our numerical experiments: synthetic test tensors and realistic applications such as the compression of facial image samples in the Olivetti database and word counts in the Enron email dataset.

Posted Content
TL;DR: Numerical experiments with transport through heterogeneous media and the Burgers' equation show orders of magnitude speedups of the proposed nonlinear reduced models based on transported subspaces compared to traditional linear reduced models and full models.
Abstract: This work presents a method for constructing online-efficient reduced models of large-scale systems governed by parametrized nonlinear scalar conservation laws. The solution manifolds induced by transport-dominated problems such as hyperbolic conservation laws typically exhibit nonlinear structures, which means that traditional model reduction methods based on linear approximations are inefficient when applied to these problems. In contrast, the approach introduced in this work derives reduced approximations that are nonlinear by explicitly composing global transport dynamics with locally linear approximations of the solution manifolds. A time-stepping scheme evolves the nonlinear reduced models by transporting local approximation spaces along the characteristic curves of the governing equations. The proposed computational procedure allows an offline/online decomposition and is online-efficient in the sense that the complexity of accurately time-stepping the nonlinear reduced model is independent of that of the full model. Numerical experiments with transport through heterogeneous media and the Burgers' equation show orders of magnitude speedups of the proposed nonlinear reduced models based on transported subspaces compared to traditional linear reduced models and full models.

Posted Content
TL;DR: Novel generalisations to entropy inequalities, multiple constraints, and kinetic energy preservation for the Euler equations are developed and tested in numerical experiments, resulting in entropy stable/dissipative grid transfers.
Abstract: For the general class of residual distribution (RD) schemes, including many finite element (such as continuous/discontinuous Galerkin) and flux reconstruction methods, an approach to construct entropy conservative/ dissipative semidiscretizations by adding suitable correction terms has been proposed by Abgrall (J.~Comp.~Phys. 372: pp. 640--666, 2018). In this work, the correction terms are characterized as solutions of certain optimization problems and are adapted to the SBP-SAT framework, focusing on discontinuous Galerkin methods. Novel generalizations to entropy inequalities, multiple constraints, and kinetic energy preservation for the Euler equations are developed and tested in numerical experiments. For all of these optimization problems, explicit solutions are provided. Additionally, the correction approach is applied for the first time to obtain a fully discrete entropy conservative/dissipative RD scheme. Here, the application of the deferred correction (DeC) method for the time integration is essential. This paper can be seen as describing a systematic method to construct structure preserving discretization, at least for the considered example.

Posted Content
TL;DR: It turns out that this simple method based on least squares regression can compete with the quasi-Monte Carlo methods in the literature which are based on lattices and digital nets.
Abstract: We consider a least squares regression algorithm for the recovery of complex-valued functions belonging to a reproducing kernel Hilbert space $H(K)$ from random data measuring the error in $L_2(D,\varrho_D)$. We prove worst-case recovery guarantees and improve on the recent new upper bounds for sampling numbers in Krieg, M. Ullrich by explicitly controlling all the involved constants with respect to the underlying spatial dimension $d$. This leads to new preasymptotic recovery bounds with high probability for the error of Hyperbolic Fourier Regression for multivariate functions. In addition, we analyze the algorithm Hyperbolic Wavelet Regression also based on least-squares to recover non-periodic functions from random samples. As a further application we reconsider the analysis of a cubature method based on plain random points with optimal weights introduced by Oettershagen. We confirm a conjecture (which was based on various numerical experiments by Oettershagen) and give improved near-optimal worst-case error bounds with high probability. It turns out that this simple method can compete with the quasi-Monte Carlo methods in the literature which are based on lattices and digital nets. Last but not least, we contribute new preasymptotic bounds to the problem of the recovery of individual functions from $n$ samples which has been already considered by Smale, Zhou; Bohn, Griebel; Cohen, Davenport, Leviatan; Chkifa, Migliorati, Nobile, Tempone; Cohen, Migliorati; Krieg and several others.

Posted Content
TL;DR: Reduced order modelling for problems for which the resulting reduced basis spaces show a slow decay of the Kolmogorov $n$-width, or, in practical calculations, its computational surrogate given by the magnitude of the eigenvalues returned by a proper orthogonal decomposition on the solution manifold is focused on.
Abstract: In this work we focus on reduced order modelling for problems for which the resulting reduced basis spaces show a slow decay of the Kolmogorov $n$-width, or, in practical calculations, its computational surrogate given by the magnitude of the eigenvalues returned by a proper orthogonal decomposition on the solution manifold. In particular, we employ an additional preprocessing during the offline phase of the reduced basis method, in order to obtain smaller reduced basis spaces. Such preprocessing is based on the composition of the snapshots with a transport map, that is a family of smooth and invertible mappings that map the physical domain of the problem into itself. Two test cases are considered: a fluid moving in a domain with deforming walls, and a fluid past a rotating cylinder. Comparison between the results of the novel offline stage and the standard one is presented.

Posted Content
TL;DR: A rigorous error analysis is performed and better convergence rates at low regularity are established than known for classical schemes in the literature so far for this new filtered low-regularity Fourier integrator for the cubic nonlinear Schrödinger equation.
Abstract: We present a new filtered low-regularity Fourier integrator for the cubic nonlinear Schrodinger equation based on recent time discretization and filtering techniques. For this new scheme, we perform a rigorous error analysis and establish better convergence rates at low regularity than known for classical schemes in the literature so far. In our error estimates, we combine the better local error properties of the new scheme with a stability analysis based on general discrete Strichartz-type estimates. The latter allow us to handle a much rougher class of solutions as the error analysis can be carried out directly at the level of $L^2$ compared to classical results \black in dimension $d$, \black which are limited to higher-order (sufficiently smooth) Sobolev spaces $H^s$ with $s>d/2$. In particular, we are able to establish a global error estimate in $L^2$ for $H^1$ solutions which is roughly of order $\tau^{ {1\over 2} + { 5-d \over 12} }$ in dimension $d \leq 3$ ($\tau$ denoting the time discretization parameter). This breaks the "natural order barrier" of $\tau^{1/2}$ for $H^1$ solutions which holds for classical numerical schemes (even in combination with suitable filter functions).

Journal ArticleDOI
TL;DR: In this article, the inner product norm preserving relaxation Runge-Kutta method is extended to general convex quantities and conservation, dissipation, or other solution properties with respect to any convex functional are enforced by the addition of a relaxation parameter.
Abstract: The framework of inner product norm preserving relaxation Runge-Kutta methods (David I. Ketcheson, \emph{Relaxation Runge-Kutta Methods: Conservation and Stability for Inner-Product Norms}, SIAM Journal on Numerical Analysis, 2019) is extended to general convex quantities. Conservation, dissipation, or other solution properties with respect to any convex functional are enforced by the addition of a {\em relaxation parameter} that multiplies the Runge-Kutta update at each step. Moreover, other desirable stability (such as strong stability preservation) and efficiency (such as low storage requirements) properties are preserved. The technique can be applied to both explicit and implicit Runge-Kutta methods and requires only a small modification to existing implementations. The computational cost at each step is the solution of one additional scalar algebraic equation for which a good initial guess is available. The effectiveness of this approach is proved analytically and demonstrated in several numerical examples, including applications to high-order entropy-conservative and entropy-stable semi-discretizations on unstructured grids for the compressible Euler and Navier-Stokes equations.