Two numerical algorithms are proposed for computing an interval matrix containing the matrix gamma function. In 2014, the author presented algorithms for enclosing all the eigenvalues and basis of invariant subspaces of $A \in \mathbb{C}^{n \times n}$. As byproducts of these algorithms, we can obtain interval matrices containing small matrices whose spectrums are included in that of $A$. In this paper, we interpret the interval matrices containing the basis and small matrices as a result of verified block diagonalization (VBD), and establish a new framework for enclosing matrix functions using the VBD. To achieve enclosure for the gamma function of the small matrices, we derive computable perturbation bounds. We can apply these bounds if input matrices satisfy conditions. We incorporate matrix argument reductions (ARs) to force the input matrices to satisfy the conditions, and develop theories for accelerating the ARs. The first algorithm uses the VBD based on a numerical spectral decomposition, and involves only cubic complexity under an assumption. The second algorithm adopts the VBD based on a numerical Jordan decomposition, and is applicable even for defective matrices. Numerical results show efficiency and robustness of the algorithms.

A popular method for finding the projection onto the intersection of two closed convex subsets in Hilbert space is Dykstra's algorithm. In this paper, we provide sufficient conditions for Dykstra's algorithm to converge rapidly, in finitely many steps. We also analyze the behaviour of Dykstra's algorithm applied to a line and a square. This case study reveals stark similarities to the method of alternating projections. Moreover, we show that Dykstra's algorithm may stall for an arbitrarily long time. Finally, we present some open problems.

N-matrices are real $n\times n$ matrices all of whose principal minors are negative. We provide (i) an $O(2^n)$ test to detect whether or not a given matrix is an N-matrix, and (ii) a characterization of N-matrices, leading to the recursive construction of every N-matrix.

This article devises a new primal-dual weak Galerkin finite element method for the convection-diffusion equation. Optimal order error estimates are established for the primal-dual weak Galerkin approximations in various discrete norms and the standard $L^2$ norms. A series of numerical experiments are conducted and reported to verify the theoretical findings.

We present the analytical formulation and the finite element solution of a fractional-order nonlocal continuum model of a Euler-Bernoulli beam. Employing consistent definitions for the fractional-order kinematic relations, the governing equations and the associated boundary conditions are derived based on variational principles. Remarkably, the fractional-order nonlocal model gives rise to a self-adjoint and positive-definite system accepting a unique solution. Further, owing to the difficulty in obtaining analytical solutions to this boundary value problem, a finite element model for the fractional-order governing equations is presented. Following a thorough validation with benchmark problems, the fractional finite element model (f-FEM) is used to study the nonlocal response of a Euler-Bernoulli beam subjected to various loading and boundary conditions. The fractional-order positive definite system will be used here to address some paradoxical results obtained for nonlocal beams through classical integral approaches to nonlocal elasticity. Although presented in the context of a 1D Euler-Bernoulli beam, the f-FEM formulation is very general and could be extended to the solution of any general fractional-order boundary value problem.

Numerous infinite dimensional dynamical systems arising in different fields have been shown to exhibit a gradient flow structure in the Wasserstein space. We construct Two Point Flux Approximation Finite Volume schemes to discretize these problems preserving their variational structure and obtaining second order accuracy in space. The choice of the discrete solver plays an important role in designing these schemes for robustness purposes. We present two applications to test the scheme and show its order of convergence.

An energy-based model of the ferroelectric polarization process is presented in the current contribution. In an energy-based setting, dielectric displacement and strain (or displacement) are the primary independent unknowns. As an internal variable, the remanent polarization vector is chosen. The model is then governed by two constitutive functions: the free energy function and the dissipation function. Choices for both functions are given. As the dissipation function for rate-independent response is non-differentiable, it is proposed to regularize the problem. Then, a variational equation can be posed, which is subsequently discretized using conforming finite elements for each quantity. We point out which kind of continuity is needed for each field (displacement, dielectric displacement and remanent polarization) is necessary to obtain a conforming method, and provide corresponding finite elements. The elements are chosen such that Gauss' law of zero charges is satisfied exactly. The discretized variational equations are solved for all unknowns at once in a single Newton iteration. We present numerical examples gained in the open source software package Netgen/NGSolve.

Big data analysis has become a crucial part of new emerging technologies such as Internet of thing (IOT), cyber-physical analysis, deep learning, anomaly detection etc. Among many other techniques, dimensionality reduction plays a key role in such analyses and facilitate the procedure of feature selection and feature extraction. Randomized algorithms are efficient tools for handling big data tensors. They accelerate decomposing large-scale data tensors by reducing the computational complexity of deterministic algorithms and also reducing the communication among different levels of memory hierarchy which is a main bottleneck in modern computing environments and architectures. In this paper, we review recent advances in randomization for computation of Tucker decomposition and Higher Order SVD (HOSVD). We discuss both random projection and sampling approaches and also single-pass and multi-pass randomized algorithms and how they can be utilized in computation of Tucker decomposition and HOSVD. Simulations on real data including weight tensors of fully connected layers of pretrained VGG-16 and VGG-19 deep neural networks and also CIFAR-10 and CIFAR-100 datasets are provided to compare performance of some of the presented algorithms.

In our work, we consider the linear least squares problem for $m\times n$-systems of linear equations $Ax = b$, $m\geq n$, such that the matrix $A$ and right-hand side vector $b$ can vary within an interval $m\times n$-matrix and an interval $m$-vector respectively. We have to compute, as sharp as possible, an interval enclosure of the set of all least squares solutions to $Ax = b$ when $A$ and $b$ independently vary within their interval bounds. Our article is devoted to the development of the so-called PPS-methods (based on Partitioning of the Parameter Set) to solve the above problem. We reduce the normal equation system, associated with the linear lest squares problem, to a special extended matrix form and produce a symmetric interval system of linear equations that is equivalent to the original interval least squares problem. To solve such symmetric system, we propose a new construction of PPS-methods, called ILSQ-PPS, which estimates the enclosure of the solution set with practical efficiency. To demonstrate the capabilities of the ILSQ-PPS method, we present a number of numerical tests and compare their results with those obtained by other methods.

We consider a Lagrangian system $L(q,\dot q) = \sum_{l=1}^{N}L^{\{l\}}(q,\dot q)$, where the $q$-variable is treated by a Generalized Additive Runge--Kutta (GARK) method. Applying the technique of discrete variations, we show how to construct symplectic schemes. Assuming the diagonal methods for the GARK given, we present some techinques for constructing the transition matrices. We address the problem of the order of the methods and discuss some semi-separable and separable problems, showing some interesting constructions of methods with non-square coefficient matrices.

Derksen proved that the spectral norm is multiplicative with respect to vertical tensor products (also known as tensor Kronecker products). We will use this result to show that the nuclear norm and other norms of interest are also multiplicative with respect to vertical tensor products. Finally, we prove an analogue of Strassen's direct sum conjecture for the nuclear norm.

A simple variant of the BDDC preconditioner in which constraints are imposed on a selected set of subobjects (subdomain subedges, subfaces and vertices between pairs of subedges) is presented. We are able to show that the condition number of the preconditioner is bounded by $C \big(1+\log (L/h)\big)^2$, where $C$ is a constant, and $h$ and $L$ are the characteristic sizes of the mesh and the subobjects, respectively. As $L$ can be chosen almost freely, the condition number can theoretically be as small as $O(1)$. We will discuss the pros and cons of the preconditioner and its application to heterogeneous problems. Numerical results on supercomputers are provided.

A slowly-varying or thin-layer multiscale assumption empowers macroscale understanding of many physical scenarios from dispersion in pipes and rivers, including beams, shells, and the modulation of nonlinear waves, to homogenisation of micro-structures. Here we begin a new exploration of the scenario where the given physics has non-local microscale interactions. We rigorously analyse the dynamics of a basic example of shear dispersion. Near each cross-section, the dynamics is expressed in the local moments of the microscale non-local effects. Centre manifold theory then supports the local modelling of the system's dynamics with coupling to neighbouring cross-sections as a non-autonomous forcing. The union over all cross-sections then provides powerful new support for the existence and emergence of a macroscale model advection-diffusion PDE global in the large, finite-sized, domain. The approach quantifies the accuracy of macroscale advection-diffusion approximations, and has the potential to open previously intractable multiscale issues to new insights.

Automatic differentiation (AD) is an ensemble of techniques that allow to evaluate accurate numerical derivatives of a mathematical function expressed in a computer programming language. In this paper we use AD for stating and solving solid mechanics problems. Given a finite element discretization of the domain, we evaluate the free energy of the solid as the integral of its strain energy density, and we make use of AD for directly obtaining the residual force vector and the tangent stiffness matrix of the problem, as the gradient and the Hessian of the free energy respectively. The result is a remarkable simplification in the statement and the solution of complex problems involving non trivial constraints systems and both geometrical and material non linearities. Together with the continuum mechanics theoretical basis, and with a description of the specific AD technique adopted, the paper illustrates the solution of a number of solid mechanics problems, with the aim of presenting a convenient numerical implementation approach, made easily available by recent programming languages, to the solid mechanics community.

A thick-restart Lanczos type algorithm is proposed for Hermitian $J$-symmetric matrices. Since Hermitian $J$-symmetric matrices possess doubly degenerate spectra or doubly multiple eigenvalues with a simple relation between the degenerate eigenvectors, we can improve the convergence of the Lanczos algorithm by restricting the search space of the Krylov subspace to that spanned by one of each pair of the degenerate eigenvector pairs. We show that the Lanczos iteration is compatible with the $J$-symmetry, so that the subspace can be split into two subspaces that are orthogonal to each other. The proposed algorithm searches for eigenvectors in one of the two subspaces without the multiplicity. The other eigenvectors paired to them can be easily reconstructed with the simple relation from the $J$-symmetry. We test our algorithm on a matrix originating from a quantum field theory.

We study the convergence rates of the semi-discrete (SD) method originally proposed in Halidias (2012), Semi-discrete approximations for stochastic differential equations and applications, International Journal of Computer Mathematics, 89(6). The SD numerical method was originally designed mainly to reproduce qualitative properties of nonlinear stochastic differential equations (SDEs). The strong convergence property of the SD method has been proved, but except for certain classes of SDEs, the order of the method was not studied. We study the order of L2-convergence and show that it can be arbitrarily close to 1/2.

Deep learning (DL) is transforming whole industries as complicated decision-making processes are being automated by Deep Neural Networks (DNNs) trained on real-world data. Driven in part by a rapidly-expanding literature on DNN approximation theory showing that DNNs can approximate a rich variety of functions, these tools are increasingly being considered for problems in scientific computing. Yet, unlike more traditional algorithms in this field, relatively little is known about DNNs from the principles of numerical analysis, namely, stability, accuracy, computational efficiency and sample complexity. In this paper we introduce a computational framework for examining DNNs in practice, and use it to study their empirical performance with regard to these issues. We examine the performance of DNNs of different widths and depths on a variety of test functions in various dimensions, including smooth and piecewise smooth functions. We also compare DL against best-in-class methods for smooth function approximation based on compressed sensing. Our main conclusion is that there is a crucial gap between the approximation theory of DNNs and their practical performance, with trained DNNs performing relatively poorly on functions for which there are strong approximation results (e.g. smooth functions), yet performing well in comparison to best-in-class methods for other functions. Finally, we present a novel practical existence theorem, which asserts the existence of a DNN architecture and training procedure which offers the same performance as current best-in-class schemes. This result indicates the potential for practical DNN approximation, and the need for future research into practical architecture design and training strategies.

Finite differences, finite elements, and their generalizations are widely used for solving partial differential equations, and their high-order variants have respective advantages and disadvantages. Traditionally, these methods are treated as different (strong vs. weak) formulations and are analyzed using different techniques (Fourier analysis or Green's functions vs. functional analysis), except for some special cases on regular grids. Recently, the authors introduced a hybrid method, called Adaptive Extended Stencil FEM or AES-FEM (Int. J. Num. Meth. Engrg., 2016, DOI:10.1002/nme.5246), which combines features of generalized finite differences and Lagrange finite elements to achieve second-order accuracy over unstructured meshes. However, its analysis was incomplete due to the lack of existing mathematical theory that unifies the formulations and analysis of these different methods. In this work, we introduce the framework of generalized weighted residuals to unify the formulation of finite differences, finite elements, and AES-FEM. In addition, we propose a unified analysis of the well-posedness, convergence, and mesh-quality dependency of these different methods. We also report numerical results with AES-FEM to verify our analysis. We show that AES-FEM improves the accuracy of generalized finite differences while reducing the mesh-quality dependency and simplifying the implementation of high-order finite elements.

A novel probabilistic numerical method for quantifying the uncertainty induced by the time integration of ordinary differential equations (ODEs) is introduced. Departing from the classical strategy to randomize ODE solvers by adding a random forcing term, we show that a probability measure over the numerical solution of ODEs can be obtained by introducing suitable random time-steps in a classical time integrator. This intrinsic randomization allows for the conservation of geometric properties of the underlying deterministic integrator such as mass conservation, symplecticity or conservation of first integrals. Weak and mean-square convergence analysis are derived. We also analyse the convergence of the Monte Carlo estimator for the proposed random time step method and show that the measure obtained with repeated sampling converges in the mean-square sense independently of the number of samples. Numerical examples including chaotic Hamiltonian systems, chemical reactions and Bayesian inferential problems illustrate the accuracy, robustness and versatility of our probabilistic numerical method.

In this paper we consider a class of convex conic programming. In particular, we propose an inexact augmented Lagrangian (I-AL) method for solving this problem, in which the augmented Lagrangian subproblems are solved approximately by a variant of Nesterov's optimal first-order method. We show that the total number of first-order iterations of the proposed I-AL method for computing an $\epsilon$-KKT solution is at most $\mathcal{O}(\epsilon^{-7/4})$. We also propose a modified I-AL method and show that it has an improved iteration-complexity $\mathcal{O}(\epsilon^{-1}\log\epsilon^{-1})$, which is so far the lowest complexity bound among all first-order I-AL type of methods for computing an $\epsilon$-KKT solution. Our complexity analysis of the I-AL methods is mainly based on an analysis on inexact proximal point algorithm (PPA) and the link between the I-AL methods and inexact PPA. It is substantially different from the existing complexity analyses of the first-order I-AL methods in the literature, which typically regard the I-AL methods as an inexact dual gradient method. Compared to the mostly related I-AL methods \cite{Lan16}, our modified I-AL method is more practically efficient and also applicable to a broader class of problems.

The paper is concerned with guaranteed a posteriori error estimates for a class of evolutionary problems related to poroelastic media governed by the quasi-static linear Biot equations. The system is decoupled employing the fixed-stress split scheme, which leads to a semi-discrete system solved iteratively. The error bounds are derived by combining a posteriori estimates for contractive mappings with those of the functional type for elliptic partial differential equations. The estimates are applicable for any approximation in the admissible functional space and are independent of the discretization method. They are fully computable, do not contain mesh dependent constants, and provide reliable global estimates of the error measured in the energy norm. Moreover, they suggest efficient error indicators for the distribution of local errors, which can be used in adaptive procedures.

We propose kernel-based collocation methods for numerical solutions to Heath-Jarrow-Morton models with Musiela parametrization. The methods can be seen as the Euler-Maruyama approximation of some finite dimensional stochastic differential equations, and allow us to compute the derivative prices by the usual Monte Carlo methods. We derive a bound on the rate of convergence under some decay condition on the inverse of the interpolation matrix and some regularity conditions on the volatility functionals.

We utilize generalized moving least squares (GMLS) to develop meshfree techniques for discretizing hydrodynamic flow problems on manifolds. We use exterior calculus to formulate incompressible hydrodynamic equations in the Stokesian regime and handle the divergence-free constraints via a generalized vector potential. This provides less coordinate-centric descriptions and enables the development of efficient numerical methods and splitting schemes for the fourth-order governing equations in terms of a system of second-order elliptic operators. Using a Hodge decomposition, we develop methods for manifolds having spherical topology. We show the methods exhibit high-order convergence rates for solving hydrodynamic flows on curved surfaces. The methods also provide general high-order approximations for the metric, curvature, and other geometric quantities of the manifold and associated exterior calculus operators. The approaches also can be utilized to develop high-order solvers for other scalar-valued and vector-valued problems on manifolds.

In this work we introduce and study novel Quasi Newton minimization methods based on a Hessian approximation Broyden Class-\textit{type} updating scheme, where a suitable matrix $\tilde{B}_k$ is updated instead of the current Hessian approximation $B_k$. We identify conditions which imply the convergence of the algorithm and, if exact line search is chosen, its quadratic termination. By a remarkable connection between the projection operation and Krylov spaces, such conditions can be ensured using low complexity matrices $\tilde{B}_k$ obtained projecting $B_k$ onto algebras of matrices diagonalized by products of two or three Householder matrices adaptively chosen step by step. Extended experimental tests show that the introduction of the adaptive criterion, which theoretically guarantees the convergence, considerably improves the robustness of the minimization schemes when compared with a non-adaptive choice; moreover, they show that the proposed methods could be particularly suitable to solve large scale problems where $L$-$BFGS$ performs poorly.

With a high probability the Sarlos algorithm of 2006 outputs a nearly optimal least squares solution of a highly overdetermined linear system of equations. We propose a simple variation that does the same for a random input with a high probability (and therefore for a large input class) rather than for any input but unlike Sarlos's algorithm performs at sublinear cost or, as we say, superfast, that is, by using much fewer memory cells and arithmetic operations than an input matrix has entries. Our extensive tests are in good accordance with this result.

The cavitation of solid elastic spheres is a classical problem of continuum mechanics. Here, we study this problem within the context of "stochastic elasticity" where the constitutive parameters are characterised by probability density functions. We consider homogeneous spheres of stochastic neo-Hookean material, composites with two concentric stochastic neo-Hookean phases, and inhomogeneous spheres of locally neo-Hookean material with a radially varying parameter. In all cases, we show that the material at the centre determines the critical load at which a spherical cavity forms there. However, while under dead-load traction, a supercritical bifurcation, with stable cavitation, is obtained in a static sphere of stochastic neo-Hookean material, for the composite and radially inhomogeneous spheres, a subcritical bifurcation, with snap cavitation, is also possible. For the dynamic spheres, oscillatory motions are produced under suitable dead-load traction, such that a cavity forms and expands to a maximum radius, then collapses again to zero periodically, but not under impulse traction. Under a surface impulse, a subcritical bifurcation is found in a static sphere of stochastic neo-Hookean material and also in an inhomogeneous sphere, whereas in composite spheres, supercritical bifurcations can occur as well. Given the non-deterministic material parameters, the results can be characterised in terms of probability distributions.

Algebraic multigrid (AMG) is a widely used scalable solver and preconditioner for large-scale linear systems resulting from the discretization of a wide class of elliptic PDEs. While AMG has optimal computational complexity, the cost of communication has become a significant bottleneck that limits its scalability as processor counts continue to grow on modern machines. This paper examines the design, implementation, and parallel performance of a novel algorithm, Algebraic Multigrid Domain Decomposition (AMG-DD), designed specifically to limit communication. The goal of AMG-DD is to provide a low-communication alternative to standard AMG V-cycles by trading some additional computational overhead for a significant reduction in communication cost. Numerical results show that AMG-DD achieves superior accuracy per communication cost compared to AMG, and speedup over AMG is demonstrated on a large GPU cluster.

In this work, we apply the (1/G')-expansion method to produce the novel soliton solution of the Gilson-Pickering equation. This method is fundamental on homogeneous balance procedure that gives the order of the estimating polynomial-type solution. Also it is based on the appreciate wave transform to reduce the governing equation. The solutions that we obtain are include of hyperbolic, complex and rational functions solutions. Finally, the results are graphically discussed.

We consider rank-1 lattices for integration and reconstruction of functions with series expansion supported on a finite index set. We explore the connection between the periodic Fourier space and the non-periodic cosine space and Chebyshev space, via tent transform and then cosine transform, to transfer known results from the periodic setting into new insights for the non-periodic settings. Fast discrete cosine transform can be applied for the reconstruction phase. To reduce the size of the auxiliary index set in the associated component-by-component (CBC) construction for the lattice generating vectors, we work with a bi-orthonormal set of basis functions, leading to three methods for function reconstruction in the non-periodic settings. We provide new theory and efficient algorithmic strategies for the CBC construction. We also interpret our results in the context of general function approximation and discrete least-squares approximation.

First-principles models of complex dynamic phenomena often have many degrees of freedom, only a small fraction of which may be scientifically relevant or observable. Reduced models distill such phenomena to their essence by modeling only relevant variables, thus decreasing computational cost and clarifying dynamical mechanisms. Here, we consider data-driven model reduction for nonlinear dynamical systems without sharp scale separation. Motivated by a discrete-time version of the Mori-Zwanzig projection operator formalism and the Wiener filter, we propose a simple and flexible mathematical formulation based on Wiener projection, which decomposes a nonlinear dynamical system into a component predictable by past values of relevant variables and its orthogonal complement. Wiener projection is equally applicable to deterministic chaotic dynamics and randomly-forced systems, and provides a natural starting point for systematic approximations. In particular, we use it to derive NARMAX models from an underlying dynamical system, thereby clarifying the scope of these widely-used tools in time series analysis. We illustrate its versatility on the Kuramoto-Sivashinsky model of spatiotemporal chaos and a stochastic Burgers equation.

Multigrid methods are popular solution algorithms for many discretized PDEs, either as standalone iterative solvers or as preconditioners, due to their high efficiency. However, the choice and optimization of multigrid components such as relaxation schemes and grid-transfer operators is crucial to the design of optimally efficient algorithms. It is well--known that local Fourier analysis (LFA) is a useful tool to predict and analyze the performance of these components. In this paper, we develop a local Fourier analysis of monolithic multigrid methods based on additive Vanka relaxation schemes for mixed finite-element discretizations of the Stokes equations. The analysis offers insight into the choice of "patches" for the Vanka relaxation, revealing that smaller patches offer more effective convergence per floating point operation. Parameters that minimize the two-grid convergence factor are proposed and numerical experiments are presented to validate the LFA predictions.

We provide a new analysis of local SGD, removing unnecessary assumptions and elaborating on the difference between two data regimes: identical and heterogeneous. In both cases, we improve the existing theory and provide values of the optimal stepsize and optimal number of local iterations. Our bounds are based on a new notion of variance that is specific to local SGD methods with different data. The tightness of our results is guaranteed by recovering known statements when we plug $H=1$, where $H$ is the number of local steps. The empirical evidence further validates the severe impact of data heterogeneity on the performance of local SGD.

We consider a function of the numerical solution of an initial value problem, its Hessian matrix with respect to the initial data, and the computation of a Hessian-vector multiplication. A simple way of approximating the Hessian-vector multiplication is to integrate the so-called second-order adjoint system numerically. However, the error in the approximation could be significant unless the numerical integration is sufficiently accurate. This paper presents a novel algorithm that computes the intended Hessian-vector multiplication exactly. For this aim, we give a new concise derivation of the second-order adjoint system and show that the intended multiplication can be computed exactly by applying a particular numerical method to the second-order adjoint system. In the discussion, symplectic partitioned Runge--Kutta methods play an important role.

We present a quantum eigenstate filtering algorithm based on quantum signal processing (QSP) and minimax polynomials. The algorithm allows us to efficiently prepare a target eigenstate of a given Hamiltonian, if we have access to an initial state with non-trivial overlap with the target eigenstate and have a reasonable lower bound for the spectral gap. We apply this algorithm to the quantum linear system problem (QLSP), and present two algorithms based on quantum adiabatic computing (AQC) and Quantum Zeno effect respectively. Both algorithms prepare the final solution as a pure state, and achieves the near optimal $\mathcal{\widetilde{O}}(d\kappa\log(1/\epsilon))$ query complexity for a $d$-sparse matrix, where $\kappa$ is the condition number, and $\epsilon$ is the desired precision. Neither algorithm uses phase estimation or amplitude amplification.

Linear multistep methods (LMMs) are popular time discretization techniques for the numerical solution of differential equations. Traditionally they are applied to solve for the state given the dynamics (the forward problem), but here we consider their application for learning the dynamics given the state (the inverse problem). This repurposing of LMMs is largely motivated by growing interest in data-driven modeling of dynamics, but the behavior and analysis of LMMs for discovery turn out to be significantly different from the well-known, existing theory for the forward problem. Assuming the highly idealized setting of being given the exact state, we establish for the first time a rigorous framework based on refined notions of consistency and stability to yield convergence using LMMs for discovery. When applying these concepts to three popular M-step LMMs, the Adams-Bashforth, Adams-Moulton, and Backwards Differentiation Formula schemes, the new theory suggests that Adams-Bashforth for $1 \leq M \leq 6$, Adams-Moulton for $M=0$ and $M=1$, and Backwards Differentiation Formula for $M \in \mathbb{N}$ are convergent, and, otherwise, the methods are not convergent in general. In addition, we provide numerical experiments to both motivate and substantiate our theoretical analysis.

Boundary element methods (BEM) reduce a partial differential equation in a domain to an integral equation on the domain's boundary. They are particularly attractive for solving problems on unbounded domains, but handling the dense matrices corresponding to the integral operators requires efficient algorithms. This article describes two approaches that allow us to solve boundary element equations on surface meshes consisting of several million of triangles while preserving the optimal convergence rates of the Galerkin discretization.

In many applications, one wants to model physical systems consisting of two different physical processes in two different domains that are coupled across a common interface. A crucial challenge is then that the solutions of the two different domains often depend critically on the interaction at the interface and therefore the problem cannot be easily decoupled into its subproblems. Here, we present a framework for finding robust preconditioners for a fairly general class of such problems by exploiting operators representing fractional and weighted Laplacians at the interface. Furthermore, we show feasibility of the framework for two common multiphysics problems; namely the Darcy-Stokes problem and a fluid--structure interaction problem. Numerical experiments that demonstrate the effectiveness of the approach are included.

The coupled Darcy-Stokes problem is widely used for modeling fluid transport in physical systems consisting of a porous part and a free part. In this work we consider preconditioners for monolitic solution algorithms of the coupled Darcy-Stokes problem, where the Darcy problem is in primal form. We employ the operator preconditioning framework and utilize a fractional solver at the interface between the problems to obtain order optimal schemes that are robust with respect to the material parameters, i.e. the permeability, viscosity and Beavers-Joseph-Saffman parameter. Our approach is similar to our earlier work, but since the Darcy problem is in primal form, the mass conservation at the interface introduces some challenges. These challenges will be specifically addressed in this paper. Numerical experiments illustrating the performance are provided. The preconditioner is posed in non-standard Sobolev spaces which may be perceived as an obstacle for its use in applications. However, we detail the implementational aspects and show that the preconditioner is quite feasible to realize in practice.

The study of greedy approximation in the context of convex optimization is becoming a promising research direction as greedy algorithms are actively being employed to construct sparse minimizers for convex functions with respect to given sets of elements. In this paper we propose a unified way of analyzing a certain kind of greedy-type algorithms for the minimization of convex functions on Banach spaces. Specifically, we define the class of Weak Biorthogonal Greedy Algorithms for convex optimization that contains a wide range of greedy algorithms. We analyze the introduced class of algorithms and establish the properties of convergence, rate of convergence, and numerical stability, which is understood in the sense that the steps of the algorithm are allowed to be performed not precisely but with controlled computational inaccuracies. We show that the following well-known algorithms for convex optimization --- the Weak Chebyshev Greedy Algorithm (co) and the Weak Greedy Algorithm with Free Relaxation (co) --- belong to this class, and introduce a new algorithm --- the Rescaled Weak Relaxed Greedy Algorithm (co). Presented numerical experiments demonstrate the practical performance of the aforementioned greedy algorithms in the setting of convex minimization as compared to optimization with regularization, which is the conventional approach of constructing sparse minimizers.

We consider stochastic differential equations driven by a general L\'evy processes (SDEs) with infinite activity and the related, via the Feynman-Kac formula, Dirichlet problem for parabolic integro-differential equation (PIDE). We approximate the solution of PIDE using a numerical method for the SDEs. The method is based on three ingredients: (i) we approximate small jumps by a diffusion; (ii) we use restricted jump-adaptive time-stepping; and (iii) between the jumps we exploit a weak Euler approximation. We prove weak convergence of the considered algorithm and present an in-depth analysis of how its error and computational cost depend on the jump activity level. Results of some numerical experiments, including pricing of barrier basket currency options, are presented.

This paper aims at an accurate and efficient computation of effective quantities, e.g., the homogenized coefficients for approximating the solutions to partial differential equations with oscillatory coefficients. Typical multiscale methods are based on a micro-macro coupling, where the macro model describes the coarse scale behaviour, and the micro model is solved only locally to upscale the effective quantities, which are missing in the macro model. The fact that the micro problems are solved over small domains within the entire macroscopic domain, implies imposing artificial boundary conditions on the boundary of the microscopic domains. A naive treatment of these artificial boundary conditions leads to a first order error in $\varepsilon/\delta$, where $\varepsilon < \delta$ represents the characteristic length of the small scale oscillations and $\delta^d$ is the size of micro domain. This error dominates all other errors originating from the discretization of the macro and the micro problems, and its reduction is a main issue in today's engineering multiscale computations. The objective of the present work is to analyze a parabolic approach, first announced in [A. Abdulle, D. Arjmand, E. Paganoni, C. R. Acad. Sci. Paris, Ser. I, 2019], for computing the homogenized coefficients with arbitrarily high convergence rates in $\varepsilon/\delta$. The analysis covers the setting of periodic micro structure, and numerical simulations are provided to verify the theoretical findings for more general settings, e.g. random stationary micro structures.

Avikainen provided a sharp upper bound of the difference $\mathbb{E}[|g(X)-g(\widehat{X})|^{q}]$ by the moments of $|X-\widehat{X}|$ for any one-dimensional random variables $X$ with bounded density and $\widehat{X}$, and function of bounded variation $g$. In this article, we generalize this estimate to any one-dimensional random variable $X$ with H\"older continuous distribution function. As applications, we provide the rate of convergence for numerical schemes for solutions of one-dimensional stochastic differential equations (SDEs) driven by Brownian motion and symmetric $\alpha$-stable with $\alpha \in (1,2)$, fractional Brownian motion with drift and Hurst parameter $H \in (0,1/2)$, and stochastic heat equations (SHEs) with Dirichlet boundary conditions driven by space--time white noise, with irregular coefficients. We also consider a numerical scheme for maximum and integral type functionals of SDEs driven by Brownian motion with irregular coefficients and payoffs which are related to multilevel Monte Carlo method.

This work concerns the numerical approximation of the multicomponent compressible Euler system for a mixture of immiscible fluids in multiple space dimensions and its contribution is twofold. We first derive an entropy stable, positive and accurate three-point finite volume scheme using relaxation-based approximate Riemann solvers from Bouchut [Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources, Frontiers in Mathematics, Birkhauser, 2004] and Coquel and Perthame [SIAM J. Numer. Anal., 35 (1998)]. Then, we extend these results to the high-order discontinuous Galerkin spectral element method (DGSEM) based on collocation of quadrature and interpolation points [Kopriva and Gassner, J. Sci. Comput., 44 (2010)]. The method relies on the framework introduced by Fisher and Carpenter [J. Comput. Phys., 252 (2013)] and Gassner [SIAM J. Sci. Comput., 35 (2013),] where we replace the physical fluxes by entropy conservative numerical fluxes [Tadmor, Math. Comput., 49 (1987)] in the integral over discretization cells, while entropy stable numerical fluxes are used at cell interfaces. Time discretization is performed with a strong-stability preserving Runge-Kutta scheme. We design two-point numerical fluxes satisfying the Tadmor's entropy conservation condition and use the numerical flux from the three-point scheme as entropy stable flux. We derive conditions on the numerical parameters to guaranty a semi-discrete entropy inequality and positivity of the fully discrete DGSEM scheme at any approximation order. The scheme is also accurate in the sense that the solution at interpolation points is exact for stationary contact waves. Numerical experiments in one and two space dimensions on flows with discontinuous solutions support the conclusions of our analysis and highlight stability, robustness and high resolution of the scheme.

In the present paper we study quasi-Monte Carlo rules for approximating integrals over the $d$-dimensional unit cube for functions from weighted Sobolev spaces of regularity one. While the properties of these rules are well understood for anchored Sobolev spaces, this is not the case for the ANOVA spaces, which are another very important type of reference spaces for quasi-Monte Carlo rules. Using a direct approach we provide a formula for the worst case error of quasi-Monte Carlo rules for functions from weighted ANOVA spaces. As a consequence we bound the worst case error from above in terms of weighted discrepancy of the employed integration nodes. On the other hand we also obtain a general lower bound in terms of the number $n$ of used integration nodes. For the one-dimensional case our results lead to the optimal integration rule and also in the two-dimensional case we provide rules yielding optimal convergence rates.

Let $f(x)$, $x\in\mathbb R^2$, be a piecewise smooth function with a jump discontinuity across a smooth surface $\mathcal S$. Let $f_{\Lambda\epsilon}$ denote the Lambda tomography (LT) reconstruction of $f$ from its discrete Radon data $\hat f(\alpha_k,p_j)$. The sampling rate along each variable is $\sim\epsilon$. First, we compute the limit $f_0(\check x)=\lim_{\epsilon\to0}\epsilon f_{\Lambda\epsilon}(x_0+\epsilon\check x)$ for a generic $x_0\in\mathcal S$. Once the limiting function $f_0(\check x)$ is known (which we call the discrete transition behavior, or DTB for short), the resolution of reconstruction can be easily found. Next, we show that straight segments of $\mathcal S$ lead to non-local artifacts in $f_{\Lambda\epsilon}$, and that these artifacts are of the same strength as the useful singularities of $f_{\Lambda\epsilon}$. We also show that $f_{\Lambda\epsilon}(x)$ does not converge to its continuous analogue $f_\Lambda=(-\Delta)^{1/2}f$ as $\epsilon\to0$ even if $x\not\in\mathcal S$. Results of numerical experiments presented in the paper confirm these conclusions. We also consider a class of Fourier integral operators $\mathcal{B}$ with the same canonical relation as the classical Radon transform adjoint, and conormal distributions $g\in\mathcal{E}'(Z_n)$, $Z_n:=S^{n-1}\times\mathbb R$, and obtain easy to use formulas for the DTB when $\mathcal{B} g$ is computed from discrete data $g(\alpha_{\vec k},p_j)$. Exact and LT reconstructions are particlular cases of this more general theory.

The absolute value equations (AVE) $Ax - |x| - b = 0$ is of interest of the optimization community. Recently, the SOR-like iteration method has been developed (Ke and Ma [{\em Appl. Math. Comput.}, 311:195--202, 2017]) and shown to be efficient for numerically solving the AVE with $\nu=\|A^{-1}\|_2<1$ (Ke and Ma [{\em Appl. Math. Comput.}, 311:195--202, 2017]; Guo, Wu and Li [{\em Appl. Math. Lett.}, 97:107--113, 2019]). Since the SOR-like iteration method is one-parameter-dependent, it is an important problem to determine the optimal iteration parameter. In this paper, we revisit the convergence conditions of the SOR-like iteration method proposed by Ke and Ma ([{\em Appl. Math. Comput.}, 311:195--202, 2017]). Furthermore, we explore the optimal parameter which minimizes $\|T(\omega)\|_2$ and the approximate optimal parameter which minimizes $\eta=\max\{|1-\omega|,\nu\omega^2\}$. The optimal and approximate optimal parameters are iteration-independent. Numerical results demonstrate that the SOR-like iteration method with the optimal parameter is superior to that with the approximate optimal parameter proposed by Guo, Wu and Li ([{\em Appl. Math. Lett.}, 97:107--113, 2019]).

Sparse iterative solvers preconditioned with the algebraic multigrid has been devised as an optimal technology to speed up the response of large linear systems. In this work, this technique was introduced into the framework of the dual delineation approach. This involves a single groundwater flow solution and a scalar advective transport solution with different right-hand sides. The method was compared with traditional preconditioned iterative methods and a direct sparse solver on several two- and three-dimensional benchmarks spanning homogeneous and heterogeneous formations. The algebraic multigrid preconditioning enabled speedups lying between one and two orders of magnitude for the groundwater flow problems. However, the sparse direct solver was the fastest for the pure advective transport processes such as the forward travel time simulations. This leads to conclude that the best sparse solver for the general advection-dispersion transport equation is P\'eclet number dependent. When equipped with the best solvers, the dual delineation technique was run on a took only a few seconds to solve multi-million grid cell problems paving the way for comprehensive sensitivity analysis. The paper gives practical hints on the strategies and conditions under which the algebraic multigrid preconditioning for nonlinear and/or transient problems would remain competitive.

A hybrid network--deterministic model for simulation of multiresistant pathogen spread in a healthcare system is presented. The model accounts for two paths of pathogen transmission between the healthcare facilities: inter-hospital patient transfers (direct transfers) and readmission of colonized patients (indirect transfers). In the latter case, the patients may be readmitted to the same facility or to a different one. Intra-hospital pathogen transmission is governed by a SIS model expressed by a system of ordinary differential equations. Using a network model created for a Lower Saxony region (Germany), we showed that the proposed model reproduces the basic properties of healthcare-associated pathogen spread. Moreover, it shows the important contribution of the readmission of colonized patients on the prevalence of individual hospitals as well as of overall healthcare system: it can increase the overall prevalence by the factor of 4 as compared to inter-hospital transfers only. The final prevalence in individual healthcare facilities was shown to depend on average length of stay by a non-linear concave function. Finally, we demonstrated that the network parameters of the model may be derived from administrative admission/discharge records. In particular, they are sufficient to obtain inter-hospital transfer probabilities, and to express the patients' transfer as a Markov process.

This paper is concerned with a numerical solution of the acoustic scattering by a bounded impenetrable obstacle in three dimensions. The obstacle scattering problem is formulated as a boundary value problem in a bounded domain by using a Dirichlet-to-Neumann (DtN) operator. An a posteriori error estimate is derived for the finite element method with the truncated DtN operator. The a posteriori error estimate consists of the finite element approximation error and the truncation error of the DtN operator, where the latter is shown to decay exponentially with respect to the truncation parameter. Based on the a posteriori error estimate, an adaptive finite element method is developed for the obstacle scattering problem. The truncation parameter is determined by the truncation error of the DtN operator and the mesh elements for local refinement are marked through the finite element approximation error. Numerical experiments are presented to demonstrate the effectiveness of the proposed method.

This paper is the first one of two serial articles, whose goal is to prove convergence of HX Preconditioner (proposed by Hiptmair and Xu 2007) for Maxwell's equations with jump coefficients. In this paper we establish various extensions of the regular Helmholtz decomposition for edge finite element functions defined in three dimensional domains. The functions defined by the regular Helmholtz decompositions can preserve the zero tangential complement on faces and edges of polyhedral domains and some non-Lipchitz domains, and possess stability estimates with only a $logarithm$ factor. These regular Helmholtz decompositions will be used to prove convergence of the HX preconditioner for Maxwell's equations with jump coefficients in another paper (Hu 2017).

Motivated by fractional derivative models in viscoelasticity, a class of semilinear stochastic Volterra integro-differential equations, and their deterministic counterparts, are considered. A generalized exponential Euler method, named here as the Mittag-Leffler Euler integrator, is used for the temporal discretization, while the spatial discretization is performed by the spectral Galerkin method. The temporal rate of strong convergence is found to be (almost) twice compared to when the backward Euler method is used together with a convolution quadrature for time discretization. Numerical experiments that validate the theory are presented.

We apply a space adaptive interior penalty discontinuous Galerkin method for solving advective Allen-Cahn equation with expanding and contracting velocity fields. The advective Allen-Cahn equation is first discretized in time and the resulting semi-linear elliptic PDE is solved by an adaptive algorithm using a residual-based a posteriori error estimator. The a posteriori error estimator contains additional terms due to the non-divergence-free velocity field. Numerical examples demonstrate the effectiveness and accuracy of the adaptive approach by resolving the sharp layers accurately.

In this paper we are concerned with numerical methods for nonhomogeneous Helmholtz equations in inhomogeneous media. We design a least squares method for discretization of the considered Helmholtz equations. In this method, an auxiliary unknown is introduced on the common interface of any two neighboring elements and a quadratic subject functional is defined by the jumps of the traces of the solutions of local Helmholtz equations across all the common interfaces, where the local Helmholtz equations are defined on elements and are imposed Robin-type boundary conditions given by the auxiliary unknowns. A minimization problem with the subject functional is proposed to determine the auxiliary unknowns. The resulting discrete system of the auxiliary unknowns is Hermitian positive definite and so it can be solved by the PCG method. Under some assumptions we show that the generated approximate solutions possess almost the optimal error estimates with little "wave number pollution". Moreover, we construct a substructuring preconditioner for the discrete system of the auxiliary unknowns. Numerical experiments show that the proposed methods are very effective for the tested Helmholtz equations with large wave numbers.

Motivated by numerical methods for solving parametric partial differential equations, this paper studies the approximation of multivariate analytic functions by algebraic polynomials. We introduce various anisotropic model classes based on Taylor expansions, and study their approximation by finite dimensional polynomial spaces $\cal{P}_{\Lambda}$ described by lower sets $\Lambda$. Given a budget $n$ for the dimension of $\cal{P}_{\Lambda}$, we prove that certain lower sets $\Lambda_n$, with cardinality $n$, provide a certifiable approximation error that is in a certain sense optimal, and that these lower sets have a simple definition in terms of simplices. Our main goal is to obtain approximation results when the number of variables $d$ is large and even infinite, and so we concentrate almost exclusively on the case $d=\infty$. We also emphasize obtaining results which hold for the full range $n\ge 1$, rather than asymptotic results that only hold for $n$ sufficiently large. In applications, one typically wants $n$ small to comply with computational budgets.

We consider the numerical approximation of the spectral fractional diffusion problem based on the so called Balakrishnan representation. The latter consists of an improper integral approximated via quadratures. At each quadrature point, a reaction-diffusion problem must be approximated and is the method bottle neck. In this work, we propose to reduce the computational cost using a reduced basis strategy allowing for a fast evaluation of the reaction-diffusion problems. The reduced basis does not depend on the fractional power $s$ for $0

Sparse spectral methods for solving partial differential equations have been derived in recent years using hierarchies of classical orthogonal polynomials on intervals, disks, and triangles. In this work we extend this methodology to a hierarchy of non-classical orthogonal polynomials on disk slices (e.g. a half-disk) and trapeziums. This builds on the observation that sparsity is guaranteed due to the boundary being defined by an algebraic curve, and that the entries of partial differential operators can be determined using formulae in terms of (non-classical) univariate orthogonal polynomials. We apply the framework to solving the Poisson, variable coefficient Helmholtz, and Biharmonic equations.

When solving linear systems arising from PDE discretizations, iterative methods (such as Conjugate Gradient, GMRES, or MINRES) are often the only practical choice. To converge in a small number of iterations, however, they have to be coupled with an efficient preconditioner. The efficiency of the preconditioner depends largely on its accuracy on the eigenvectors corresponding to small eigenvalues, and unfortunately, black-box methods typically cannot guarantee sufficient accuracy on these eigenvectors. Thus, constructing the preconditioner becomes a problem-dependent task. However, for a large class of problems, including many elliptic equations, the eigenvectors corresponding to small eigenvalues are smooth functions of the PDE grid. In this paper, we describe a hierarchical approximate factorization approach which focuses on improving accuracy on the smooth eigenvectors. The improved accuracy is achieved by preserving the action of the factorized matrix on piecewise polynomial functions of the grid. Based on the factorization, we propose a family of sparse preconditioners with $O(n)$ or $O(n \log{n})$ construction complexities. Our methods exhibit the optimal $O(n)$ solution times in benchmarks run on large elliptic problems of different types, arising for example in flow or mechanical simulations. In the case of the linear elasticity equation the preconditioners are exact on the near-kernel rigid body modes.

In this paper, authors focus effort on improving the conventional discrete velocity method (DVM) into a multiscale scheme in finite volume framework for gas flow in all flow regimes. Unlike the typical multiscale kinetic methods unified gas-kinetic scheme (UGKS) and discrete unified gas-kinetic scheme (DUGKS), which concentrate on the evolution of the distribution function at the cell interface, in the present scheme the flux for macroscopic variables is split into the equilibrium part and the nonequilibrium part, and the nonequilibrium flux is calculated by integrating the discrete distribution function at the cell center, which overcomes the excess numerical dissipation of the conventional DVM in the continuum flow regime. Afterwards, the macroscopic variables are finally updated by simply integrating the discrete distribution function at the cell center, or by a blend of the increments based on the macroscopic and the microscopic systems, and the multiscale property is achieved. Several test cases, involving unsteady, steady, high speed, low speed gas flows in all flow regimes, have been performed, demonstrating the good performance of the multiscale DVM from free molecule to continuum Navier-Stokes solutions and the multiscale property of the scheme is proved.

Methods for quantifying the effects of uncertainties in hyperbolic problems can be divided into intrusive and non-intrusive techniques. Non-intrusive methods allow the usage of a given deterministic solver in a black-box manner, while being embarrassingly parallel. However, avoiding intrusive modifications of a given solver takes away the ability to use several inherently intrusive numerical acceleration tools. Moreover, intrusive methods are expected to reach a given accuracy with a smaller number of unknowns compared to non-intrusive techniques. This effect is amplified in settings with high dimensional uncertainty. A downside of intrusive methods is however the need to guarantee hyperbolicity of the resulting moment system. In contrast to stochastic-Galerkin (SG), the Intrusive Polynomial Moment (IPM) method is able to maintain hyperbolicity at the cost of solving an optimization problem in every spatial cell and every time step. In this work, we propose several acceleration techniques for intrusive methods and study their advantages and shortcomings compared to the non-intrusive Stochastic Collocation method. When solving steady problems with IPM, the numerical costs arising from repeatedly solving the IPM optimization problem can be reduced by using concepts from PDE-constrained optimization. Additionally, we propose an adaptive implementation and efficient parallelization strategy of the IPM method. The effectiveness of the proposed adaptations is demonstrated for multi-dimensional uncertainties in fluid dynamics applications, resulting in the observation of requiring a smaller number of unknowns to achieve a given accuracy when using intrusive methods. Furthermore, using the proposed acceleration techniques, our implementation reaches a given accuracy faster than Stochastic Collocation.

Many problems in materials science and biology involve particles interacting with strong, short-ranged bonds, that can break and form on experimental timescales. Treating such bonds as constraints can significantly speed up sampling their equilibrium distribution, and there are several methods to sample subject to fixed constraints. We introduce a Monte Carlo method to handle the case when constraints can break and form. Abstractly, the method samples a probability distribution on a stratification: a collection of manifolds of different dimensions, where the lower-dimensional manifolds lie on the boundaries of the higher-dimensional manifolds. We show several applications of the method in polymer physics, self-assembly of colloids, and volume calculation in high dimensions.