Abstract: The Lambert W-function is the solution to the transcendental equation W(x)e W(x) = x. It has two real branches, one of which, for x ∈ [−1/e, ∞], is usually denoted as the principal branch. On this branch, the function grows from − 1 to infinity, logarithmically at large x. The present work is devoted to the construction of accurate approximations for the principal branch of the W-function. In particular, a simple, global analytic approximation is derived that covers the whole branch with a maximum relative error smaller than 5 × 10−3. Starting from it, machine precision accuracy is reached everywhere with only three steps of a quadratically convergent iterative scheme, here examined for the first time, which is more efficient than standard Newton’s iteration at large x. Analytic bounds for W are also constructed, for x > e, which are much tighter than those currently available. It is noted that the exponential of the upper bounding function yields an upper bound for the prime counting function π(n) that is better than the well-known Chebyshev’s estimates at large n. Finally, the construction of accurate approximations to W based on Chebyshev spectral theory is discussed; the difficulties involved are highlighted, and methods to overcome them are presented.

Abstract: A two-point boundary value problem is considered on the interval [0, 1], where the leading term in the differential operator is a Riemann-Liouville fractional derivative of order 2 ź ź with 0 <ź < 1. It is shown that any solution of such a problem can be expressed in terms of solutions to two associated weakly singular Volterra integral equations of the second kind. As a consequence, existence and uniqueness of a solution to the boundary value problem are proved, the structure of this solution is elucidated, and sharp bounds on its derivatives (in terms of the parameter ź) are derived. These results show that in general the first-order derivative of the solution will blow up at x = 0, so accurate numerical solution of this class of problems is not straightforward. The reformulation of the boundary problem in terms of Volterra integral equations enables the construction of two distinct numerical methods for its solution, both based on piecewise-polynomial collocation. Convergence rates for these methods are proved and numerical results are presented to demonstrate their performance.

Abstract: In boundary integral methods it is often necessary to evaluate layer potentials on or close to the boundary, where the underlying integral is difficult to evaluate numerically. Quadrature by expansion (QBX) is a new method for dealing with such integrals, and it is based on forming a local expansion of the layer potential close to the boundary. In doing so, one introduces a new quadrature error due to nearly singular integration in the evaluation of expansion coefficients. Using a method based on contour integration and calculus of residues, the quadrature error of nearly singular integrals can be accurately estimated. This makes it possible to derive accurate estimates for the quadrature errors related to QBX, when applied to layer potentials in two and three dimensions. As examples we derive estimates for the Laplace and Helmholtz single layer potentials. These results can be used for parameter selection in practical applications.

Abstract: We use asymptotically optimal adaptive numerical methods (here specifically a wavelet scheme) for snapshot computations within the offline phase of the Reduced Basis Method (RBM). The resulting discretizations for each snapshot (i.e., parameter-dependent) do not permit the standard RB `truth space', but allow for error estimation of the RB approximation with respect to the exact solution of the considered parameterized partial differential equation. The residual-based a posteriori error estimators are computed by an adaptive dual wavelet expansion, which allows us to compute a surrogate of the dual norm of the residual. The resulting adaptive RBM is analyzed. We show the convergence of the resulting adaptive greedy method. Numerical experiments for stationary and instationary problems underline the potential of this approach.

Abstract: In this article, we apply the newly introduced numerical method which is a combination of Sumudu transforms and Homotopy analysis method for the solution of time fractional third order dispersive type PDE equations. It is also discussed generalized algorithm, absolute convergence and analytic result of the finite number of independent variables including time variable.

Abstract: Tensor-compressed numerical solution of elliptic multiscale-diffusion and high frequency scattering problems is considered. For either problem class, solutions exhibit multiple length scales governed by the corresponding scale parameter: the scale of oscillations of the diffusion coefficient or smallest wavelength, respectively. As is well-known, this imposes a scale-resolution requirement on the number of degrees of freedom required to accurately represent the solutions in standard finite-element (FE) discretizations. Low-order FE methods are by now generally perceived unsuitable for high-frequency coefficients in diffusion problems and high wavenumbers in scattering problems. Accordingly, special techniques have been proposed instead (such as numerical homogenization, heterogeneous multiscale method, oversampling, etc.) which require, in some form, a-priori information on the microstructure of the solution. We analyze the approximation properties of tensor-formatted, conforming first-order FE methods for scale resolution in multiscale problems without a-priori information. The FE methods are based on the dynamic extraction of principal components from stiffness matrices, load and solution vectors by the quantized tensor train (QTT) decomposition. For prototypical model problems, we prove that this approach, by means of the QTT reparametrization of the FE space, allows to identify effective degrees of freedom to replace the degrees of freedom of a uniform "virtual" (i.e. never directly accessed) mesh, whose number may be prohibitively large to realize computationally. Precisely, solutions of model elliptic homogenization and high-frequency acoustic scattering problems are proved to admit QTT-structured approximations whose number of effective degrees of freedom required to reach a prescribed approximation error scales polylogarithmically with respect to the reciprocal of the target Sobolev-norm accuracy ź with only a mild dependence on the scale parameter. No a-priori information on the nature of the problems and intrinsic length scales of the solution is required in the numerical realization of the presently proposed QTT-structured approach. Although only univariate model multiscale problems are analyzed in the present paper, QTT structured algorithms are applicable also in several variables. Detailed numerical experiments confirm the theoretical bounds. As a corollary of our analysis, we prove that for the mentioned model problems, the Kolmogorov n-widths of solution sets are exponentially small for analytic data, independently of the problems' scale parameters. That implies, in particular, the exponential convergence of reduced basis techniques which is scale-robust, i.e., independent of the scale parameter in the problem.

Abstract: Recently, there has been considerable work on analyzing learning algorithms with pairwise loss functions in the batch setting. There is relatively little theoretical work on analyzing their online algorithms, despite of their popularity in practice due to the scalability to big data. In this paper, we consider online learning algorithms with pairwise loss functions based on regularization schemes in reproducing kernel Hilbert spaces. In particular, we establish the convergence of the last iterate of the online algorithm under a very weak assumption on the step sizes and derive satisfactory convergence rates for polynomially decaying step sizes. Our technique uses Rademacher complexities which handle function classes associated with pairwise loss functions. Since pairwise learning involves pairs of examples, which are no longer i.i.d., standard techniques do not directly apply to such pairwise learning algorithms. Hence, our results are a non-trivial extension of those in the setting of univariate loss functions to the pairwise setting.

Abstract: In this paper we generalize and improve a recently developed domain decomposition preconditioner for the iterative solution of discretized Helmholtz equations. We introduce an improved method for transmission at the internal boundaries using perfectly matched layers. Simultaneous forward and backward sweeps are introduced, thereby improving the possibilities for parallellization. Finally, the method is combined with an outer two-grid iteration. The method is studied theoretically and with numerical examples. It is shown that the modifications lead to substantial decreases in computation time and memory use, so that computation times become comparable to that of the fastests methods currently in the literature for problems with up to 108 degrees of freedom.

Abstract: In this paper, we analyze the supercloseness property of the streamline diffusion finite element method (SDFEM) on Shishkin triangular meshes, which is different from one in the case of rectangular meshes. The analysis depends on integral inequalities for the parts related to the diffusion in the bilinear form. Moreover, our result allows the construction of a simple postprocessing that yields a more accurate solution. Finally, numerical experiments support these theoretical results.

Abstract: A rotation---minimizing frame (f1,f2,f3) on a space curve r(ź) defines an orthonormal basis for ź3$\mathbb {R}^{3}$ in which f1=rź/|rź|$\mathbf {f}_{1}=\mathbf {r}^{\prime }/|\mathbf {r}^{\prime }|$ is the curve tangent, and the normal---plane vectors f2, f3 exhibit no instantaneous rotation about f1. Polynomial curves that admit rational rotation---minimizing frames (or RRMF curves) form a subset of the Pythagorean---hodograph (PH) curves, specified by integrating the form rź(ź)=A(ź)iAź(ź)$\mathbf {r}^{\prime }(\xi )=\mathcal {A}(\xi )\,\mathbf{i} \,\mathcal {A}^{*}(\xi )$ for some quaternion polynomial A(ź)$\mathcal {A}(\xi )$. By introducing the notion of the rotation indicatrix and the core of the quaternion polynomial A(ź)$\mathcal {A}(\xi )$, a comprehensive characterization of the complete space of RRMF curves is developed, that subsumes all previously known special cases. This novel characterization helps clarify the structure of the complete space of RRMF curves, distinguishes the spatial RRMF curves from trivial (planar) cases, and paves the way toward new construction algorithms.

Abstract: The classic eigenvalue problem of the Laplace operator inside a variety of polygons is numerically solved by using a method nearly identical to that used by Fox, Henrici, and Moler in their 1967 paper. It is demonstrated that such eigenvalue calculations can be extended to unprecedented precision, often to well over a hundred digits, or even thousands of digits. To work well, geometric symmetry must be exploited. The de-symmetrized fundamental domains (usually triangular) considered here have at most one non-analytic vertex. Dirichlet, Neumann, and periodic-type edge conditions are independently imposed on each symmetry-reduced polygon edge. The method of particular solutions is used whereby an eigenfunction is expanded in an N-term Fourier-Bessel series about the non-analytic vertex and made to match at a set of N points on the boundary. Under the right conditions, the so-called point-matching determinant has roots that approximate eigenvalues. A key observation is that by increasing the number of terms in the expansion, the approximate eigenvalue may be made to alternate above and below, while approaching what is presumed to be the exact eigenvalue. This alternation effectively provides a new method to bound eigenvalues, by inspection. Specific examples include Dirichlet and Neumann eigenvalues within polygons with re-entrant angles (L-shape, cut-square, 5-point star) and the regular polygons. Thousand-digit results are reported for the lowest Dirichlet eigenvalues of the L-shape, and regular pentagon and hexagon.

Abstract: We consider the problem of spatiotemporal sampling in a discrete infinite dimensional spatially invariant evolutionary process x(n) = Anx to recover an unknown convolution operator A given by a filter a?l1(?)$a \in \ell ^{1}(\mathbb {Z})$ and an unknown initial state x modeled as a vector in l2(?)$\ell ^{2}(\mathbb {Z})$. Traditionally, under appropriate hypotheses, any x can be recovered from its samples on ?$\mathbb {Z}$ and A can be recovered by the classical techniques of deconvolution. In this paper, we will exploit the spatiotemporal correlation and propose a new sampling scheme to recover A and x that allows us to sample the evolving states x,Ax,? ,AN?1x on a sub-lattice of ?$\mathbb {Z}$, and thus achieve a spatiotemporal trade off. The spatiotemporal trade off is motivated by several industrial applications (Lu and Vetterli, 2249---2252, 2009). Specifically, we show that {x(m?),Ax(m?),?,AN?1x(m?):N?2m}$\{x(m\mathbb {Z}), Ax(m\mathbb {Z}), \cdots , A^{N-1}x(m\mathbb {Z}): N \geq 2m\}$ contains enough information to recover a typical "low pass filter" a and x almost surely, thus generalizing the idea of the finite dimensional case in Aldroubi and Krishtal, arXiv:1412.1538 (2014). In particular, we provide an algorithm based on a generalized Prony method for the case when both a and x are of finite impulse response and an upper bound of their support is known. We also perform a perturbation analysis based on the spectral properties of the operator A and initial state x, and verify the results by several numerical experiments. Finally, we provide several other numerical techniques to stabilize the proposed method, with some examples to demonstrate the improvement.

Abstract: Linear and semilinear second-order parabolic equations are considered. For these equations, we give a posteriori error estimates in the maximum norm that improve upon recent results in the literature. In particular it is shown that logarithmic dependence on the time step size can be eliminated. Semidiscrete and fully discrete versions of the backward Euler and of the Crank-Nicolson methods are considered. For their full discretizations, we use elliptic reconstructions that are, respectively, piecewise-constant and piecewise-linear in time. Certain bounds for the Green’s function of the parabolic operator are also employed.

Abstract: In this paper, we construct a second order algorithm based on the spectral deferred correction method for the time-dependent magnetohydrodynamics flows at a low magnetic Reynolds number. We present a complete theoretical analysis to prove that this algorithm is unconditionally stable, consistent and second order accuracy. Finally, two numerical examples are given to illustrate the convergence and effectiveness of our algorithm.

Abstract: The goal of this paper is to introduce a simple finite element method to solve the Stokes equations. This method is in primal velocity-pressure formulation and is so simple such that both velocity and pressure are approximated by piecewise constant functions. Implementation issues as well as error analysis are investigated. A basis for a divergence free subspace of the velocity field is constructed so that the original saddle point problem can be reduced to a symmetric and positive definite system with much fewer unknowns. The numerical experiments indicate that the method is accurate.

Abstract: Generalized shift-invariant (GSI) systems, originally introduced by Hernandez et al. and Ron and Shen, provide a common frame work for analysis of Gabor systems, wavelet systems, wave packet systems, and other types of structured function systems. In this paper we analyze three important aspects of such systems. First, in contrast to the known cases of Gabor frames and wavelet frames, we show that for a GSI system forming a frame, the Calderon sum is not necessarily bounded by the lower frame bound. We identify a technical condition implying that the Calderon sum is bounded by the lower frame bound and show that under a weak assumption the condition is equivalent with the local integrability condition introduced by Hernandez et al. Second, we provide explicit and general constructions of frames and dual pairs of frames having the GSI-structure. In particular, the setup applies to wave packet systems and in contrast to the constructions in the literature, these constructions are not based on characteristic functions in the Fourier domain. Third, our results provide insight into the local integrability condition (LIC).

Abstract: Families of parameter dependent univariate and bivariate subdivision schemes are presented in this paper. These families are new variants of the Lane-Riesenfeld algorithm. So the subdivision algorithms consist of both refining and smoothing steps. In refining step, we use the quartic B-spline based subdivision schemes. In smoothing step, we average the adjacent points. The bivariate schemes are the non-tensor product version of our univariate schemes. Moreover, for odd and even number of smoothing steps, we get the primal and dual schemes respectively. Higher regularity of the schemes can be achieved by increasing the number of smoothing steps. These schemes can be nicely generalized to contain local shape parameters that allow the user to adjust locally the shape of the limit curve/surface.

Abstract: We give an algorithm for approximating a given plane curve segment by a planar elastic curve. The method depends on an analytic representation of the space of elastic curve segments, together with a geometric method for obtaining a good initial guess for the approximating curve. A gradient-driven optimization is then used to find the approximating elastic curve.

Abstract: We propose a new family of finite volume high-accurate numerical schemes devoted to solve one-dimensional steady-state hyperbolic systems. High-accuracy (up to the sixth-order presently) is achieved thanks to polynomial reconstructions while stability is provided with an a posteriori MOOD method which control the cell polynomial degree for eliminating non-physical oscillations in the vicinity of discontinuities. Such a procedure demands the determination of a chain detector to discriminate between troubled and valid cells, a cascade of polynomial degrees to be successively tested when oscillations are detected, and a parachute scheme corresponding to the last, viscous, and robust scheme of the cascade. Experimented on linear, Burgers', and Euler equations, we demonstrate that the schemes manage to retrieve smooth solutions with optimal order of accuracy but also irregular solutions without spurious oscillations.

Abstract: In this paper we propose and analyze a new family of nonlinear subdivision schemes which can be considered non-oscillatory versions of the 6-point Deslauries-Dubuc (DD) interpolatory scheme, just as the Power p schemes are considered nonlinear non-oscillatory versions of the 4-point DD interpolatory scheme. Their design principle may be related to that of the Power p schemes and it is based on a weighted analog of the Power p mean. We prove that the new schemes reproduce exactly polynomials of degree three and stay ’close’ to the 6-point DD scheme in smooth regions. In addition, we prove that the first and second difference schemes are well defined for each member of the family, which allows us to give a simple proof of the uniform convergence of these schemes and also to study their stability as in [19, 22]. However our theoretical study of stability is not conclusive and we perform a series of numerical experiments that seem to point out that only a few members of the new family of schemes are stable. On the other hand, extensive numerical testing reveals that, for smooth data, the approximation order and the regularity of the limit function may be similar to that of the 6-point DD scheme and larger than what is obtained with the Power p schemes.

Abstract: Minimal Polynomial Extrapolation (MPE) and Reduced Rank Extrapolation (RRE) are two polynomial methods used for accelerating the convergence of sequences of vectors {xm}. They are applied successfully in conjunction with fixed-point iterative schemes in the solution of large and sparse systems of linear and nonlinear equations in different disciplines of science and engineering. Both methods produce approximations sk to the limit or antilimit of {xm} that are of the form sk=źi=0kźixi$\boldsymbol {s}_{k}={\sum }^{k}_{i=0}\gamma _{i}\boldsymbol {x}_{i}$ with źi=0kźi=1${\sum }^{k}_{i=0}\gamma _{i}=1$, for some scalars źi. The way the two methods are derived suggests that they might, somehow, be related to each other; this has not been explored so far, however. In this work, we tackle this issue and show that the vectors skMPE$\boldsymbol {s}_{k}^{\textit {{\tiny {MPE}}}}$ and skRRE$\boldsymbol {s}_{k}^{\textit {{\tiny {RRE}}}}$ produced by the two methods are related in more than one way, and independently of the way the xm are generated. One of our results states that RRE stagnates, in the sense that skRRE=skź1RRE$\boldsymbol {s}_{k}^{\textit {{\tiny {RRE}}}}=\boldsymbol {s}_{k-1}^{\textit {{\tiny {RRE}}}}$, if and only if skMPE$\boldsymbol {s}_{k}^{\textit {{\tiny {MPE}}}}$ does not exist. Another result states that, when skMPE$\boldsymbol {s}_{k}^{\textit {{\tiny {MPE}}}}$ exists, there holds μkskRRE=μkź1skź1RRE+źkskMPEwithμk=μkź1+źk,$$\mu_{k}\boldsymbol{s}_{k}^{\textit{{\tiny{RRE}}}}=\mu_{k-1}\boldsymbol{s}_{k-1}^{\textit{{\tiny{RRE}}}}+ u_{k}\boldsymbol{s}_{k}^{\textit{{\tiny{MPE}}}}\quad \text{with}\quad \mu_{k}=\mu_{k-1}+ u_{k}, $$for some positive scalars μk, μkź1, and źk that depend only on skRRE$\boldsymbol {s}_{k}^{\textit {{\tiny {RRE}}}}$, skź1RRE$\boldsymbol {s}_{k-1}^{\textit {{\tiny {RRE}}}}$, and skMPE$\boldsymbol {s}_{k}^{\textit {{\tiny {MPE}}}}$, respectively. Our results are valid when MPE and RRE are defined in any weighted inner product and the norm induced by it. They also contain as special cases the known results pertaining to the connection between the method of Arnoldi and the method of generalized minimal residuals, two important Krylov subspace methods for solving nonsingular linear systems.

Abstract: In this paper, a class of high-order central Hermite WENO (HWENO) schemes based on finite volume framework and staggered meshes is proposed for directly solving one- and two-dimensional Hamilton-Jacobi (HJ) equations. The methods involve the Lax-Wendroff type discretizations or the natural continuous extension of Runge-Kutta methods in time. This work can be regarded as an extension of central HWENO schemes for hyperbolic conservation laws (Tao et al. J. Comput. Phys. 318, 222–251, 2016) which combine the central scheme and the HWENO spatial reconstructions and therefore carry many features of both schemes. Generally, it is not straightforward to design a finite volume scheme to directly solve HJ equations and a key ingredient for directly solving such equations is the reconstruction of numerical Hamiltonians to guarantee the stability of methods. Benefited from the central strategy, our methods require no numerical Hamiltonians. Meanwhile, the zeroth-order and the first-order moments of the solution are involved in the spatial HWENO reconstructions which is more compact compared with WENO schemes. The reconstructions are implemented through a dimension-by-dimension strategy when the spatial dimension is higher than one. A collection of one- and two- dimensional numerical examples is performed to validate high resolution and robustness of the methods in approximating the solutions of HJ equations, which involve linear, nonlinear, smooth, non-smooth, convex or non-convex Hamiltonians.

Abstract: The main purpose of this paper is to investigate the piecewise spectral collocation method for system of Volterra integral equations. The provided convergence analysis shows that the presented method performs better than global spectral collocation method and piecewise polynomial collocation method. Numerical experiments are carried out to confirm these theoretical results.

Abstract: The paper provides an approach for constructing multivariate radial kernels satisfying higher-order generalized Strang-Fix conditions from a given univariate generator. There are three key features of the approach. First, the kernels are explicitly expressed only by the derivatives of the f-form of the generator without computing any Fourier transforms. Second, it includes the radial kernels with the highest-order generalized Strang-Fix conditions. Finally, it requires only computing univariate derivatives of the f-form. Therefore, the approach is simple, efficient and easy to implement. As examples, the paper constructs radial kernels from some commonly used generators, including the Gaussian functions, the inverse multiquadric functions and compactly supported positive definite functions.

Abstract: A well-balanced van Leer-type numerical scheme for the shallow water equations with variable topography is presented. The model involves a nonconservative term, which often makes standard schemes difficult to approximate solutions in certain regions. The construction of our scheme is based on exact solutions in computational form of local Riemann problems. Numerical tests are conducted, where comparisons between this van Leer-type scheme and a Godunov-type scheme are provided. Data for the tests are taken in both the subcritical region as well as supercritical region. Especially, tests for resonant cases where the exact solutions contain coinciding waves are also investigated. All numerical tests show that each of these two methods can give a good accuracy, while the van Leer -type scheme gives a better accuracy than the Godunov-type scheme. Furthermore, it is shown that the van Leer-type scheme is also well-balanced in the sense that it can capture exactly stationary contact discontinuity waves.

Abstract: A model order reduction method is developed for an operator with a non-empty null-space and applied to numerical solution of a forward multi-frequency eddy current problem using a rational interpolation of the transfer function in the complex plane. The equation is decomposed into the part in the null space of the operator, calculated exactly, and the part orthogonal to it which is approximated on a low-dimensional rational Krylov subspace. For the Maxwell's equations the null space is related to the null space of the curl. The proposed null space correction is related to divergence correction and uses the Helmholtz decomposition. In the case of the finite element discretization with the edge elements, it is accomplished by solving the Poisson equation on the nodal elements of the same grid. To construct the low-dimensional approximation we adaptively choose the interpolating frequencies, defining the rational Krylov subspace, to reduce the maximal approximation error. We prove that in the case of an adaptive choice of shifts, the matrix spanning the approximation subspace can never become rank deficient. The efficiency of the developed approach is demonstrated by applying it to the magnetotelluric problem, which is a geophysical electromagnetic remote sensing method used in mineral, geothermal, and groundwater exploration. Numerical tests show an excellent performance of the proposed methods characterized by a significant reduction of the computational time without a loss of accuracy. The null space correction regularizes the otherwise ill-posed interpolation problem.

Abstract: We consider an h-p version of the continuous Petrov-Galerkin time stepping method for Volterra integro-differential equations with proportional delays. We derive a priori error bounds in the L 2-, H 1- and L ∞-norm that are explicit in the local time steps, the local approximation orders, and the local regularity of the exact solution. Numerical experiments are presented to illustrate the theoretical results.

Abstract: This paper presents an integral formulation for Helmholtz problems with mixed boundary conditions. Unlike most integral equation techniques for mixed boundary value problems, the proposed method uses a global boundary charge density. As a result, Calderon identities can be utilized to avoid the use of hypersingular integral operators. Numerical results illustrate the performance of the proposed solution technique.

Abstract: Because physical phenomena on Earth's surface occur on many different length scales, it makes sense when seeking an efficient approximation to start with a crude global approximation, and then make a sequence of corrections on finer and finer scales. It also makes sense eventually to seek fine scale features locally, rather than globally. In the present work, we start with a global multiscale radial basis function (RBF) approximation, based on a sequence of point sets with decreasing mesh norm, and a sequence of (spherical) radial basis functions with proportionally decreasing scale centered at the points. We then prove that we can "zoom in" on a region of particular interest, by carrying out further stages of multiscale refinement on a local region. The proof combines multiscale techniques for the sphere from Le Gia, Sloan and Wendland, SIAM J. Numer. Anal. 48 (2010) and Applied Comp. Harm. Anal. 32 (2012), with those for a bounded region in ?d from Wendland, Numer. Math. 116 (2010). The zooming in process can be continued indefinitely, since the condition numbers of matrices at the different scales remain bounded. A numerical example illustrates the process.