Authors:David Brander; Jens Gravesen; Toke Bjerge Nørbjerg Pages: 25 - 43 Abstract: 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. PubDate: 2017-02-01 DOI: 10.1007/s10444-016-9474-z Issue No:Vol. 43, No. 1 (2017)

Authors:Christiaan C. Stolk Pages: 45 - 76 Abstract: 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. PubDate: 2017-02-01 DOI: 10.1007/s10444-016-9475-y Issue No:Vol. 43, No. 1 (2017)

Authors:Natalia Kopteva; Martin Stynes Pages: 77 - 99 Abstract: 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. PubDate: 2017-02-01 DOI: 10.1007/s10444-016-9476-x Issue No:Vol. 43, No. 1 (2017)

Authors:Christian Gerhards; Sergiy Pereverzyev; Pavlo Tkachenko Pages: 101 - 112 Abstract: Abstract In many geoscientific applications, multiple noisy observations of different origin need to be combined to improve the reconstruction of a common underlying quantity. This naturally leads to multi-parameter models for which adequate strategies are required to choose a set of ‘good’ parameters. In this study, we present a fairly general method for choosing such a set of parameters, provided that discrete direct, but maybe noisy, measurements of the underlying quantity are included in the observation data, and the inner product of the reconstruction space can be accurately estimated by the inner product of the discretization space. Then the proposed parameter choice method gives an accuracy that only by an absolute constant multiplier differs from the noise level and the accuracy of the best approximant in the reconstruction and in the discretization spaces. PubDate: 2017-02-01 DOI: 10.1007/s10444-016-9477-9 Issue No:Vol. 43, No. 1 (2017)

Authors:François Bachoc; Martin Ehler; Manuel Gräf Pages: 113 - 126 Abstract: Abstract Motivated by the construction of confidence intervals in statistics, we study optimal configurations of 2 d − 1 lines in real projective space ℝℙ d−1. For small d, we determine line sets that numerically minimize a wide variety of potential functions among all configurations of 2 d − 1 lines through the origin. Numerical experiments verify that our findings enable to assess efficiently the tightness of a bound arising from the statistical literature. PubDate: 2017-02-01 DOI: 10.1007/s10444-016-9478-8 Issue No:Vol. 43, No. 1 (2017)

Authors:Avram Sidi Pages: 151 - 170 Abstract: Abstract Minimal Polynomial Extrapolation (MPE) and Reduced Rank Extrapolation (RRE) are two polynomial methods used for accelerating the convergence of sequences of vectors {x m }. 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 s k to the limit or antilimit of {x m } that are of the form \(\boldsymbol {s}_{k}={\sum }^{k}_{i=0}\gamma _{i}\boldsymbol {x}_{i}\) with \({\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 \(\boldsymbol {s}_{k}^{\textit {{\tiny {MPE}}}}\) and \(\boldsymbol {s}_{k}^{\textit {{\tiny {RRE}}}}\) produced by the two methods are related in more than one way, and independently of the way the x m are generated. One of our results states that RRE stagnates, in the sense that \(\boldsymbol {s}_{k}^{\textit {{\tiny {RRE}}}}=\boldsymbol {s}_{k-1}^{\textit {{\tiny {RRE}}}}\) , if and only if \(\boldsymbol {s}_{k}^{\textit {{\tiny {MPE}}}}\) does not exist. Another result states that, when \(\boldsymbol {s}_{k}^{\textit {{\tiny {MPE}}}}\) exists, there holds $$\mu_{k}\boldsymbol{s}_{k}^{\textit{{\tiny{RRE}}}}=\mu_{k-1}\boldsymbol{s}_{k-1}^{\textit{{\tiny{RRE}}}}+ \nu_{k}\boldsymbol{s}_{k}^{\textit{{\tiny{MPE}}}}\quad \text{with}\quad \mu_{k}=\mu_{k-1}+\nu_{k}, $$ for some positive scalars μ k , μ k−1, and ν k that depend only on \(\boldsymbol {s}_{k}^{\textit {{\tiny {RRE}}}}\) , \(\boldsymbol {s}_{k-1}^{\textit {{\tiny {RRE}}}}\) , and \(\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. PubDate: 2017-02-01 DOI: 10.1007/s10444-016-9481-0 Issue No:Vol. 43, No. 1 (2017)

Authors:Michal Kordy; Elena Cherkaev; Philip Wannamaker Pages: 171 - 193 Abstract: 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. PubDate: 2017-02-01 DOI: 10.1007/s10444-016-9482-z Issue No:Vol. 43, No. 1 (2017)

Authors:Ludvig af Klinteberg; Anna-Karin Tornberg Pages: 195 - 234 Abstract: 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. PubDate: 2017-02-01 DOI: 10.1007/s10444-016-9484-x Issue No:Vol. 43, No. 1 (2017)

Authors:Ghulam Mustafa; Rabia Hameed Abstract: 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. PubDate: 2017-03-13 DOI: 10.1007/s10444-017-9519-y

Authors:Pietro Dell’Acqua Abstract: Abstract In recent years, several efforts were made in order to introduce boundary conditions for deblurring problems that allow to get accurate reconstructions. This resulted in the birth of Reflective, Anti-Reflective and Mean boundary conditions, which are all based on the idea of guaranteeing the continuity of the signal/image outside the boundary. Here we propose new boundary conditions that are obtained by suitably combining Taylor series and finite difference approximations. Moreover, we show that also Anti-Reflective and Mean boundary conditions can be attributed to the same framework. Numerical results show that, in case of low levels of noise and blurs able to perform a suitable smoothing effect on the original image (e.g. Gaussian blur), the proposed boundary conditions lead to a significant improvement of the restoration accuracy with respect to those available in the literature. PubDate: 2017-03-10 DOI: 10.1007/s10444-017-9525-0

Authors:Fu-Jun Liu; Zhong-Qing Wang; Hui-Yuan Li Abstract: Abstract A fully diagonalized spectral method using generalized Laguerre functions is proposed and analyzed for solving elliptic equations on the half line. We first define the generalized Laguerre functions which are complete and mutually orthogonal with respect to an equivalent Sobolev inner product. Then the Fourier-like Sobolev orthogonal basis functions are constructed for the diagonalized Laguerre spectral method of elliptic equations. Besides, a unified orthogonal Laguerre projection is established for various elliptic equations. On the basis of this orthogonal Laguerre projection, we obtain optimal error estimates of the fully diagonalized Laguerre spectral method for both Dirichlet and Robin boundary value problems. Finally, numerical experiments, which are in agreement with the theoretical analysis, demonstrate the effectiveness and the spectral accuracy of our diagonalized method. PubDate: 2017-03-08 DOI: 10.1007/s10444-017-9522-3

Authors:William Paulsen; Samuel Cowgill Abstract: Abstract The generalized tetration, defined by the equation F(z+1) = b F(z) in the complex plane with F(0) = 1, is considered for any b > e 1/e . By comparing other solutions to Kneser’s solution, natural conditions are found which force Kneser’s solution to be the unique solution to the equation. This answers a conjecture posed by Trappmann and Kouznetsov. Also, a new iteration method is developed which numerically approximates the function F(z) with an error of less than 10−50 for most bases b, using only 180 nodes, with each iteration gaining one or two places of accuracy. This method can be applied to other problems involving the Abel equation. PubDate: 2017-03-07 DOI: 10.1007/s10444-017-9524-1

Authors:Natalia Kopteva; Torsten Linß Abstract: 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. PubDate: 2017-03-02 DOI: 10.1007/s10444-017-9514-3

Authors:Theresa Wenger; Sina Ober-Blöbaum; Sigrid Leyendecker Abstract: Abstract In this work, variational integrators of higher order for dynamical systems with holonomic constraints are constructed and analyzed. The construction is based on approximating the configuration and the Lagrange multiplier via different polynomials. The splitting of the augmented Lagrangian in two parts enables the use of different quadrature formulas to approximate the integral of each part. Conditions are derived that ensure the linear independence of the higher order constrained discrete Euler-Lagrange equations and stiff accuracy. Time reversibility is investigated for the discrete flow on configuration level only as for the flow on configuration and momentum level. The fulfillment of the hidden constraints plays an important role for the time reversibility of the presented integrators. The order of convergence is investigated numerically. Order reduction of the momentum and the Lagrange multiplier compared to the order of the configuration occurs in general, but can be avoided by fulfilling the hidden constraints in a simple post processing step. Regarding efficiency versus accuracy a numerical analysis yields that higher orders increase the accuracy of the discrete solution substantially while the computational costs decrease. A comparison to the constrained Galerkin methods in Marsden and West (Acta Numerica 10, 357–514 2001) and the symplectic SPARK integrators of Jay (SIAM Journal on Numerical Analysis 45(5), 1814–1842 2007) reveals that the approach presented here is more general and thus allows for more flexibility in the design of the integrator. PubDate: 2017-03-02 DOI: 10.1007/s10444-017-9520-5

Authors:X. Claeys; R. Hiptmair; E. Spindler Abstract: Abstract We consider isotropic scalar diffusion boundary value problems whose diffusion coefficients are piecewise constant with respect to a partition of space into Lipschitz subdomains. We allow so-called material junctions where three or more subdomains may abut. We derive a boundary integral equation of the second kind posed on the skeleton of the subdomain partition that involves, as unknown, only one trace function at each point of each interface. We prove the well-posedness of the corresponding boundary integral equations. We also report numerical tests for Galerkin boundary element discretisations, in which the new approach proves to be highly competitive compared to the well-established first kind direct single-trace boundary integral formulation. In particular, GMRES seems to enjoy fast convergence independent of the mesh resolution for the discrete second kind BIE. PubDate: 2017-03-01 DOI: 10.1007/s10444-017-9517-0

Authors:J. A. Ezquerro; M. A. Hernández-Verón Abstract: Abstract This paper focuses on the importance of center conditions on the first derivative of the operator involved in the solution of nonlinear equations by Newton’s method when the semilocal convergence of the method is established from the technique of recurrence relations. PubDate: 2017-02-28 DOI: 10.1007/s10444-017-9518-z

Authors:Dao Huy Cuong; Mai Duc Thanh Abstract: 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. PubDate: 2017-02-23 DOI: 10.1007/s10444-017-9521-4

Authors:Rosa Donat; Sergio López-Ureña; Maria Santágueda Abstract: 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. PubDate: 2017-02-14 DOI: 10.1007/s10444-016-9509-5

Authors:Zhanjing Tao; Jianxian Qiu Abstract: 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. PubDate: 2017-02-04 DOI: 10.1007/s10444-017-9515-2

Authors:Yanzhao Cao; Ying Jiang; Yuesheng Xu Abstract: Abstract The goal of this paper is to construct an efficient numerical algorithm for computing the coefficient matrix and the right hand side of the linear system resulting from the spectral Galerkin approximation of a stochastic elliptic partial differential equation. We establish that the proposed algorithm achieves an exponential convergence with requiring only O \((n\log _{2}^{d+1}n)\) number of arithmetic operations, where n is the highest degree of the one dimensional orthogonal polynomial used in the algorithm, d+1 is the number of terms in the finite Karhunen–Loéve (K-L) expansion. Numerical experiments confirm the theoretical estimates of the proposed algorithm and demonstrate its computational efficiency. PubDate: 2017-02-04 DOI: 10.1007/s10444-017-9513-4