Abstract: Abstract In this paper, we consider a numerical approximation of the Van Roosbroeck’s drift–diffusion system given by a backward Euler in time and finite volume in space discretization, with Scharfetter–Gummel fluxes. We first propose a proof of existence of a solution to the scheme which does not require any assumption on the time step. The result relies on the application of a topological degree argument which is based on the positivity and on uniform-in-time upper bounds of the approximate densities. Secondly, we establish uniform-in-time lower bounds satisfied by the approximate densities. These uniform-in-time upper and lower bounds ensure the exponential decay of the scheme towards the thermal equilibrium as shown in Bessemoulin-Chatard (Numer Math 25(3):147–168, 2016). PubDate: 2019-04-01

Abstract: Abstract We consider the preconditioned conjugate gradient method (PCG) with optimal preconditioner in the frame of the boundary element method for elliptic first-kind integral equations. Our adaptive algorithm steers the termination of PCG as well as the local mesh-refinement. Besides convergence with optimal algebraic rates, we also prove almost optimal computational complexity. In particular, we provide an additive Schwarz preconditioner which can be computed in linear complexity and which is optimal in the sense that the condition numbers of the preconditioned systems are uniformly bounded. As model problem serves the 2D or 3D Laplace operator and the associated weakly-singular integral equation with energy space \(\widetilde{H}^{-1/2}(\Gamma )\) . The main results also hold for the hyper-singular integral equation with energy space \(H^{1/2}(\Gamma )\) . PubDate: 2019-04-01

Abstract: Abstract The objective of this article is to characterize the entropy and \(L_2\) stability of several representative discontinuous Galerkin (DG) methods for solving the compressible Euler equations. Towards this end, three DG methods are constructed: one DG method with entropy variables as its unknowns, and two DG methods with conservative variables as their unknowns. These methods are employed in order to discretize the compressible Euler equations in space. Thereafter, the resulting semi-discrete formulations are analyzed, and the entropy and \(L_2\) stability characteristics are evaluated. It is shown that the semi-discrete formulation of the DG method with entropy variables is entropy and \(L_2\) stable. Furthermore, it is shown that the semi-discrete formulations of the DG methods with conservative variables are only guaranteed to be entropy and \(L_2\) stable under the following assumptions: the entropy projection errors vanish, or the terms containing the entropy projection errors are non-positive. Thereafter, the semi-discrete formulation with entropy variables, and one of the semi-discrete formulations with conservative variables, are discretized in time with an ‘algebraically stable’ Runge–Kutta (RK) scheme. The resulting formulations are fully-discrete and can be immediately applied to practical problems. In this article, they are employed to simulate a vortex propagating for long distances. It is shown that temporal stability is maintained by the DG method with entropy variables, but the DG method with conservative variables exhibits instability. PubDate: 2019-04-01

Abstract: Abstract We consider the use of complete radiation boundary conditions for the solution of the Helmholtz equation in waveguides. A general analysis of well-posedness, convergence, and finite element approximation is given. In addition, methods for the optimization of the boundary condition parameters are considered. The theoretical results are illustrated by some simple numerical experiments. PubDate: 2019-04-01

Abstract: Abstract The biharmonic operator plays a central role in a wide array of physical models, such as elasticity theory and the streamfunction formulation of the Navier–Stokes equations. Its spectral theory has been extensively studied. In particular the one-dimensional case (over an interval) constitutes the basic model of a high order Sturm–Liouville problem. The need for corresponding numerical simulations has led to numerous works. The present paper relies on a discrete biharmonic calculus. The primary object of this calculus is a high-order compact discrete biharmonic operator (DBO). The DBO is constructed in terms of the discrete Hermitian derivative. However, the underlying reason for its accuracy remained unclear. This paper is a contribution in this direction, expounding the strong connection between cubic spline functions (on an interval) and the DBO. The first observation is that the (scaled) fourth-order distributional derivative of the cubic spline is identical to the action of the DBO on grid functions. It is shown that the kernel of the inverse of the discrete operator is (up to scaling) equal to the grid evaluation of the kernel of \(\left[ \left( \frac{d}{dx}\right) ^4\right] ^{-1}\) , and explicit expressions are presented for both kernels. As an important application, the relation between the (infinite) set of eigenvalues of the fourth-order Sturm–Liouville problem and the finite set of eigenvalues of the discrete biharmonic operator is studied. The discrete eigenvalues are proved to converge (at an “optimal” \(O(h^4)\) rate) to the continuous ones. Another consequence is the validity of a comparison principle. It is well known that there is no maximum principle for the fourth-order equation. However, a positivity result is derived, both for the continuous and the discrete biharmonic equation, showing that in both cases the kernels are order preserving. PubDate: 2019-04-01

Abstract: Abstract Fully discrete Galerkin finite element methods are studied for the equations of miscible displacement in porous media with the commonly-used Bear–Scheidegger diffusion–dispersion tensor: $$\begin{aligned} D(\mathbf{u}) = \gamma d_m I + \mathbf{u} \bigg ( \alpha _T I + (\alpha _L - \alpha _T) \frac{\mathbf{u} \otimes \mathbf{u}}{ \mathbf{u} ^2}\bigg ) \, . \end{aligned}$$ Previous works on optimal-order \(L^\infty (0,T;L^2)\) -norm error estimate required the regularity assumption \(\nabla _x\partial _tD(\mathbf{u}(x,t)) \in L^\infty (0,T;L^\infty (\Omega ))\) , while the Bear–Scheidegger diffusion–dispersion tensor is only Lipschitz continuous even for a smooth velocity field \(\mathbf{u}\) . In terms of the maximal \(L^p\) -regularity of fully discrete finite element solutions of parabolic equations, optimal error estimate in \(L^p(0,T;L^q)\) -norm and almost optimal error estimate in \(L^\infty (0,T;L^q)\) -norm are established under the assumption of \(D(\mathbf{u})\) being Lipschitz continuous with respect to \(\mathbf{u}\) . PubDate: 2019-04-01

Abstract: Abstract Numerical approximation of a stochastic partial integro-differential equation driven by a space-time white noise is studied by truncating a series representation of the noise, with finite element method for spatial discretization and convolution quadrature for time discretization. Sharp-order convergence of the numerical solutions is proved up to a logarithmic factor. Numerical examples are provided to support the theoretical analysis. PubDate: 2019-04-01

Abstract: The corrected version states the parameter range as p ∈ 2ℕ instead of p ∈ ℕ. The effect is to disallow non-smooth norms such as the ℓ1-norm for the distance measure. PubDate: 2019-03-21

Abstract: Abstract We introduce a discontinuous Galerkin method for the mixed formulation of the elasticity eigenproblem with reduced symmetry. The analysis of the resulting discrete eigenproblem does not fit in the standard spectral approximation framework since the underlying source operator is not compact and the scheme is nonconforming. We show that the proposed scheme provides a correct approximation of the spectrum and prove asymptotic error estimates for the eigenvalues and the eigenfunctions. Finally, we provide several numerical tests to illustrate the performance of the method and confirm the theoretical results. PubDate: 2019-03-20

Abstract: Abstract We consider and analyze applying a spectral inverse iteration algorithm and its subspace iteration variant for computing eigenpairs of an elliptic operator with random coefficients. With these iterative algorithms the solution is sought from a finite dimensional space formed as the tensor product of the approximation space for the underlying stochastic function space, and the approximation space for the underlying spatial function space. Sparse polynomial approximation is employed to obtain the first one, while classical finite elements are employed to obtain the latter. An error analysis is presented for the asymptotic convergence of the spectral inverse iteration to the smallest eigenvalue and the associated eigenvector of the problem. A series of detailed numerical experiments supports the conclusions of this analysis. PubDate: 2019-03-19

Abstract: Abstract We consider curvature depending variational models for image regularization, such as Euler’s elastica. These models are known to provide strong priors for the continuity of edges and hence have important applications in shape- and image processing. We consider a lifted convex representation of these models in the roto-translation space: in this space, curvature depending variational energies are represented by means of a convex functional defined on divergence free vector fields. The line energies are then easily extended to any scalar function. It yields a natural generalization of the total variation to curvature-dependent energies. As our main result, we show that the proposed convex representation is tight for characteristic functions of smooth shapes. We also discuss cases where this representation fails. For numerical solution, we propose a staggered grid discretization based on an averaged Raviart–Thomas finite elements approximation. This discretization is consistent, up to minor details, with the underlying continuous model. The resulting non-smooth convex optimization problem is solved using a first-order primal-dual algorithm. We illustrate the results of our numerical algorithm on various problems from shape- and image processing. PubDate: 2019-03-16

Abstract: Abstract Markov expanding maps, a class of simple chaotic systems, are commonly used as models for chaotic dynamics, but existing numerical methods to study long-time statistical properties such as invariant measures have a poor trade-off between computational effort and accuracy. We develop a spectral Galerkin method for these maps’ transfer operators, estimating statistical quantities using finite submatrices of the transfer operators’ infinite Fourier or Chebyshev basis coefficient matrices. Rates of convergence of these estimates are obtained via quantitative bounds on the full transfer operator matrix entries; we find the method furnishes up to exponentially accurate estimates of statistical properties in only a polynomially large computational time. To implement these results we suggest and demonstrate two algorithms: a rigorously-validated algorithm, and a fast, more convenient adaptive algorithm. Using the first algorithm we prove rigorous bounds on some exemplar quantities that are substantially more accurate than previous. We show that the adaptive algorithm can produce double floating-point accuracy estimates in a fraction of a second on a personal computer. PubDate: 2019-03-07

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

Abstract: Abstract We present natural axisymmetric variants of schemes for curvature flows introduced earlier by the present authors and analyze them in detail. Although numerical methods for geometric flows have been used frequently in axisymmetric settings, numerical analysis results so far are rare. In this paper, we present stability, equidistribution, existence and uniqueness results for the introduced approximations. Numerical computations show that these schemes are very efficient in computing numerical solutions of geometric flows as only a spatially one-dimensional problem has to be solved. The good mesh properties of the schemes also allow them to compute in very complex axisymmetric geometries. PubDate: 2019-03-01

Abstract: Abstract The discontinuous Galerkin (dG) methodology provides a hierarchy of time discretization schemes for evolutionary problems such as elastoplasticity with the Prandtl-Reuß flow rule. A dG time discretization has been proposed for a variational inequality in the context of rate-independent inelastic material behaviour in Alberty and Carstensen in (CMME 191:4949–4968, 2002) with the help of duality in convex analysis to justify certain jump terms. This paper establishes the first a priori error analysis for the dG(1) scheme with discontinuous piecewise linear polynomials in the temporal and lowest-order finite elements for the spatial discretization. Compared to a generalized mid-point rule, the dG(1) formulation distributes the action of the material law in the form of the variational inequality in time and so it introduces an error in the material law. This may result in a suboptimal convergence rate for the dG(1) scheme and this paper shows that the stress error in the \(L^\infty (L^2)\) norm is merely \(O(h+k^{3/2})\) based on a seemingly sharp error analysis. The numerical investigation for a benchmark problem with known analytic solution provides empirical evidence of a higher convergence rate of the dG(1) scheme compared to dG(0). PubDate: 2019-03-01

Abstract: Abstract This paper is concerned with the approximation of tensors using tree-based tensor formats, which are tensor networks whose graphs are dimension partition trees. We consider Hilbert tensor spaces of multivariate functions defined on a product set equipped with a probability measure. This includes the case of multidimensional arrays corresponding to finite product sets. We propose and analyse an algorithm for the construction of an approximation using only point evaluations of a multivariate function, or evaluations of some entries of a multidimensional array. The algorithm is a variant of higher-order singular value decomposition which constructs a hierarchy of subspaces associated with the different nodes of the tree and a corresponding hierarchy of interpolation operators. Optimal subspaces are estimated using empirical principal component analysis of interpolations of partial random evaluations of the function. The algorithm is able to provide an approximation in any tree-based format with either a prescribed rank or a prescribed relative error, with a number of evaluations of the order of the storage complexity of the approximation format. Under some assumptions on the estimation of principal components, we prove that the algorithm provides either a quasi-optimal approximation with a given rank, or an approximation satisfying the prescribed relative error, up to constants depending on the tree and the properties of interpolation operators. The analysis takes into account the discretization errors for the approximation of infinite-dimensional tensors. For a tensor with finite and known rank in a tree-based format, the algorithm is able to recover the tensor in a stable way using a number of evaluations equal to the storage complexity of the representation of the tensor in this format. Several numerical examples illustrate the main results and the behavior of the algorithm for the approximation of high-dimensional functions using hierarchical Tucker or tensor train tensor formats, and the approximation of univariate functions using tensorization. PubDate: 2019-03-01

Abstract: Abstract We present a new algorithm which, given a bidiagonal decomposition of a totally nonnegative matrix, computes all its eigenvalues to high relative accuracy in floating point arithmetic in \(O(n^3)\) time. It also computes exactly the Jordan blocks corresponding to zero eigenvalues in up to \(O(n^4)\) time. PubDate: 2019-03-01

Abstract: Abstract Several important problems in partial differential equations can be formulated as integral equations. Often the integral operator defines the solution of an elliptic problem with specified jump conditions at an interface. In principle the integral equation can be solved by replacing the integral operator with a finite difference calculation on a regular grid. A practical method of this type has been developed by the second author. In this paper we prove the validity of a simplified version of this method for the Dirichlet problem in a general domain in \({\mathbb {R}}^2\) or \({\mathbb {R}}^3\) . Given a boundary value, we solve for a discrete version of the density of the double layer potential using a low order interface method. It produces the Shortley–Weller solution for the unknown harmonic function with accuracy \(O(h^2)\) . We prove the unique solvability for the density, with bounds in norms based on the energy or Dirichlet norm, using techniques which mimic those of exact potentials. The analysis reveals that this crude method maintains much of the mathematical structure of the classical integral equation. Examples are included. PubDate: 2019-03-01

Abstract: Abstract We consider commutator-free exponential integrators as put forward in Alverman and Fehske (J Comput Phys 230:5930–5956, 2011). For parabolic problems, it is important for the well-definedness that such an integrator satisfies a positivity condition such that essentially it only proceeds forward in time. We prove that this requirement implies maximal convergence order of four for real coefficients, which has been conjectured earlier by other authors. PubDate: 2019-03-01

Abstract: Abstract This paper is concerned with diffusive approximations of some numerical schemes for several linear (or weakly nonlinear) kinetic models which are motivated by wide-range applications, including radiative transfer or neutron transport, run-and-tumble models of chemotaxis dynamics, and Vlasov–Fokker–Planck plasma modeling. The well-balanced method applied to such kinetic equations leads to time-marching schemes involving a “scattering S-matrix”, itself derived from a normal modes decomposition of the stationary solution. One common feature these models share is the type of diffusive approximation: their macroscopic densities solve drift-diffusion systems, for which a distinguished numerical scheme is Il’in/Scharfetter–Gummel’s “exponential fitting” discretization. We prove that these well-balanced schemes relax, within a parabolic rescaling, towards such type of discretization by means of an appropriate decomposition of the S-matrix, hence are asymptotic preserving. PubDate: 2019-03-01