Authors:Terence J. T. Norton Abstract: Abstract In choosing a numerical method for the long-time integration of reversible Hamiltonian systems one must take into consideration several key factors: order of the method, ability to preserve invariants of the system, and efficiency of the computation. In this paper, 6th-order composite symmetric general linear methods (COSY-GLMs) are constructed using a generalisation of the composition theory associated with Runge–Kutta methods (RKMs). A novel aspect of this approach involves a nonlinear transformation which is used to convert the GLM to a canonical form in which its starting and finishing methods are trivial. Numerical experiments include efficiency comparisons to symmetric diagonally-implicit RKMs, where it is shown that COSY-GLMs of the same order typically require half the number of function evaluations, as well as long-time computations of both separable and non-separable Hamiltonian systems which demonstrate the preservation properties of the new methods. PubDate: 2018-01-10 DOI: 10.1007/s10543-017-0692-7

Authors:Caterina Fenu; Lothar Reichel; Giuseppe Rodriguez; Hassane Sadok Pages: 1019 - 1039 Abstract: Abstract Tikhonov regularization is commonly used for the solution of linear discrete ill-posed problems with error-contaminated data. A regularization parameter that determines the quality of the computed solution has to be chosen. One of the most popular approaches to choosing this parameter is to minimize the Generalized Cross Validation (GCV) function. The minimum can be determined quite inexpensively when the matrix A that defines the linear discrete ill-posed problem is small enough to rapidly compute its singular value decomposition (SVD). We are interested in the solution of linear discrete ill-posed problems with a matrix A that is too large to make the computation of its complete SVD feasible, and show how upper and lower bounds for the numerator and denominator of the GCV function can be determined fairly inexpensively for large matrices A by computing only a few of the largest singular values and associated singular vectors of A. These bounds are used to determine a suitable value of the regularization parameter. Computed examples illustrate the performance of the proposed method. PubDate: 2017-12-01 DOI: 10.1007/s10543-017-0662-0 Issue No:Vol. 57, No. 4 (2017)

Authors:Davoud Mirzaei Pages: 1041 - 1063 Abstract: Abstract In this paper a direct approximation method on the sphere, constructed by generalized moving least squares, is presented and analyzed. It is motivated by numerical solution of partial differential equations on spheres and other manifolds. The new method generalizes the finite difference methods, someway, for scattered data points on each local subdomain. As an application, the Laplace–Beltrami equation is solved and the theoretical and experimental results are given. The new approach eliminates some drawbacks of the previous methods. PubDate: 2017-12-01 DOI: 10.1007/s10543-017-0659-8 Issue No:Vol. 57, No. 4 (2017)

Authors:Huiqing Xie Pages: 1065 - 1082 Abstract: Abstract A new method is proposed to compute the eigenvector derivative of a quadratic eigenvalue problem (QEP) analytically dependent on a parameter. It avoids the linearization of the QEP. The proposed method can be seen as an improved incomplete modal method. Only a few eigenvectors of the QEP are required. The contributions of other eigenvectors to the desired eigenvector derivative are obtained by an iterative scheme. From this point of view, our method also can be seen as an iterative method. The convergence properties of the proposed method are analyzed. The techniques to accelerate the proposed method are provided. A strategy is developed for simultaneously computing several eigenvector derivatives by the proposed method. Finally some numerical examples are given to demonstrate the efficiency of our method. PubDate: 2017-12-01 DOI: 10.1007/s10543-017-0680-y Issue No:Vol. 57, No. 4 (2017)

Authors:Ching-Sung Liu; Che-Rung Lee Pages: 1083 - 1108 Abstract: Abstract In this paper, we present the inexact structure preserving Arnoldi methods for structured eigenvalue problems. They are called structure preserving because the computed eigenvalues and eigenvectors can preserve the desirable properties of the structures of the original matrices, even with large errors involved in the computation of matrix-vector products. A backward error matrix is called structured if it has the same structure as the original matrix. We derive a common form for the structured backward errors that can be used for different structure preserving processes, and prove the derived form has the minimum Frobenius norm among all possible backward errors. Furthermore, we employ the derived backward errors for some specific structure preserving processes to estimate the accuracy of the solutions obtained by inexact Arnoldi methods for eigenvalue problems. We aim to give, wherever possible, formulae that are inexpensive to compute so that they can be used routinely in practice. Numerical experiments are provided to support the theoretical results. PubDate: 2017-12-01 DOI: 10.1007/s10543-017-0660-2 Issue No:Vol. 57, No. 4 (2017)

Authors:Iwona Skalna; Milan Hladík Pages: 1109 - 1136 Abstract: Abstract We propose a new approach to computing a parametric solution (the so-called p-solution) to parametric interval linear systems. Solving such system is an important part of many scientific and engineering problems involving uncertainties. The parametric solution has many useful properties. It permits to compute an outer solution, an inner estimate of the interval hull solution, and intervals containing the lower and upper bounds of the interval hull solution. It can also be employed for solving various constrained optimisation problems related to the parametric interval linear system. The proposed approach improves both outer and inner bounds for the parametric solution set. In this respect, the new approach is competitive to most of the existing methods for solving parametric interval linear systems. Improved bounds on the parametric solution set guarantees improved bounds for the solutions of related optimisation problems. PubDate: 2017-12-01 DOI: 10.1007/s10543-017-0679-4 Issue No:Vol. 57, No. 4 (2017)

Authors:Yuji Nakatsukasa Pages: 1137 - 1152 Abstract: Abstract The standard approach to computing an approximate SVD of a large-scale matrix is to project it onto lower-dimensional trial subspaces from both sides, compute the SVD of the small projected matrix, and project it back to the original space. This results in a low-rank approximate SVD to the original matrix, and we can then obtain approximate left and right singular subspaces by extracting subsets from the approximate SVD. In this work we assess the quality of the extraction process in terms of the accuracy of the approximate singular subspaces, measured by the angle between the exact and extracted subspaces (relative to the angle between the exact and trial subspaces). The main message is that the extracted approximate subspaces are optimal usually to within a modest constant. PubDate: 2017-12-01 DOI: 10.1007/s10543-017-0665-x Issue No:Vol. 57, No. 4 (2017)

Authors:Patrick Zulian; Teseo Schneider; Kai Hormann; Rolf Krause Pages: 1185 - 1203 Abstract: Abstract The discretization of the computational domain plays a central role in the finite element method. In the standard discretization the domain is triangulated with a mesh and its boundary is approximated by a polygon. The boundary approximation induces a geometry-related error which influences the accuracy of the solution. To control this geometry-related error, iso-parametric finite elements and iso-geometric analysis allow for high order approximation of smooth boundary features. We present an alternative approach which combines parametric finite elements with smooth bijective mappings leaving the choice of approximation spaces free. Our approach allows to represent arbitrarily complex geometries on coarse meshes with curved edges, regardless of the domain boundary complexity. The main idea is to use a bijective mapping for automatically warping the volume of a simple parameterization domain to the complex computational domain, thus creating a curved mesh of the latter. Numerical examples provide evidence that our method has lower approximation error for domains with complex shapes than the standard finite element method, because we are able to solve the problem directly on the exact domain without having to approximate it. PubDate: 2017-12-01 DOI: 10.1007/s10543-017-0669-6 Issue No:Vol. 57, No. 4 (2017)

Authors:Sverre Anmarkrud; Kristian Debrabant; Anne Kværnø Abstract: Abstract In this paper stochastic partitioned Runge–Kutta (SPRK) methods are considered. A general order theory for SPRK methods based on stochastic B-series and multicolored, multishaped rooted trees is developed. The theory is applied to prove the order of some known methods, and it is shown how the number of order conditions can be reduced in some special cases, especially that the conditions for preserving quadratic invariants can be used as simplifying assumptions. PubDate: 2017-12-22 DOI: 10.1007/s10543-017-0693-6

Authors:J. M. Carnicer; C. Godés Abstract: Abstract The Lebesgue constant is a measure for the stability of the Lagrange interpolation. The decomposition of the Lagrange interpolation operator in their even and odd parts with respect to the last variable can be used to find a relation between the Lebesgue constant for a space of polynomials and the corresponding Lebesgue constants for subspaces of even and odd polynomials. It is shown that such a decomposition preserves the stability properties of the Lagrange interpolation operator. We use the Lebesgue functions to provide pointwise quantitative measures of the stability properties and illustrate with examples the behaviour in simple cases. PubDate: 2017-12-01 DOI: 10.1007/s10543-017-0674-9

Authors:Jeremy Schmitt; Tatiana Shingel; Melvin Leok Abstract: Abstract In this paper, we present a variational integrator that is based on an approximation of the Euler–Lagrange boundary-value problem via Taylor’s method. This can be viewed as a special case of the shooting-based variational integrator. The Taylor variational integrator exploits the structure of the Taylor method, which results in a shooting method that is one order higher compared to other shooting methods based on a one-step method of the same order. In addition, this method can generate quadrature nodal evaluations at the cost of a polynomial evaluation, which may increase its efficiency relative to other shooting-based variational integrators. A symmetric version of the method is proposed, and numerical experiments are conducted to exhibit the efficacy and efficiency of the method. PubDate: 2017-11-24 DOI: 10.1007/s10543-017-0690-9

Authors:Marcel Klinge; Rüdiger Weiner; Helmut Podhaisky Abstract: Abstract In this paper, explicit peer methods are studied in which some of the stage values are copies of stage values from previous steps. This allows to reduce the number of function calls per step and can be interpreted as being a generalization of the first-same-as-last principle known from Runge–Kutta methods. The variable step size implementation is more complex as the nodes depend on the history of previous step size changes. Optimally zero stable explicit peer methods up to order \(p=8\) are constructed using constraint numerical optimization. In addition the constructed methods are superconvergent of order \(s+1\) for constant step sizes. The new methods show their efficiency in comparison with the Matlab codes ode23, ode45 and ode113 in numerical experiments. PubDate: 2017-11-22 DOI: 10.1007/s10543-017-0691-8

Authors:Abtin Rahimian; Alex Barnett; Denis Zorin Abstract: Abstract We introduce a quadrature scheme—QBKIX —for the ubiquitous high-order accurate evaluation of singular layer potentials associated with general elliptic PDEs, i.e., a scheme that yields high accuracy at all distances to the domain boundary as well as on the boundary itself. Relying solely on point evaluations of the underlying kernel, our scheme is essentially PDE-independent; in particular, no analytic expansion nor addition theorem is required. Moreover, it applies to boundary integrals with singular, weakly singular, and hypersingular kernels. Our work builds upon quadrature by expansion, which approximates the potential by an analytic expansion in the neighborhood of each expansion center. In contrast, we use a sum of fundamental solutions lying on a ring enclosing the neighborhood, and solve a small dense linear system for their coefficients to match the potential on a smaller concentric ring. We test the new method with Laplace, Helmholtz, Yukawa, Stokes, and Navier (elastostatic) kernels in two dimensions (2D) using adaptive, panel-based boundary quadratures on smooth and corner domains. Advantages of the algorithm include its relative simplicity of implementation, immediate extension to new kernels, dimension-independence (allowing simple generalization to 3D), and compatibility with fast algorithms such as the kernel-independent FMM. PubDate: 2017-11-06 DOI: 10.1007/s10543-017-0689-2

Authors:Jeonghun J. Lee Abstract: Abstract We propose a new finite element method for a three-field formulation of Biot’s consolidation model in poroelasticity and prove the a priori error estimates. Uniform-in-time error estimates of all the unknowns are obtained for both semidiscrete solutions and fully discrete solutions with the backward Euler time discretization. In contrast to previous results, the implicit constants in our error estimates are uniformly bounded as the Lamé coefficient indicating incompressiblity of poroelastic medium is arbitrarily large, and as the constrained specific storage coefficient is arbitrarily small. Therefore the method does not suffer from the volumetric locking of linear elasticity and provides robust error estimates without additional assumptions on the constrained specific storage coefficient. PubDate: 2017-10-22 DOI: 10.1007/s10543-017-0688-3

Authors:Jana Burkotová; Irena Rachůnková; Ewa B. Weinmüller Abstract: Abstract This paper deals with the collocation method applied to solve systems of singular linear ordinary differential equations with variable coefficient matrices and nonsmooth inhomogeneities. The classical stage convergence order is shown to hold for the piecewise polynomial collocation applied to boundary value problems with time singularities of the first kind provided that their solutions are appropriately smooth. The convergence theory is illustrated by numerical examples. PubDate: 2017-10-16 DOI: 10.1007/s10543-017-0686-5

Authors:Pavel Strachota; Michal Beneš Abstract: Abstract The Allen–Cahn equation originates in the phase field formulation of phase transition phenomena. It is a reaction-diffusion ODE with a nonlinear reaction term which allows the formation of a diffuse phase interface. We first introduce a model initial boundary-value problem for the isotropic variant of the equation. Its numerical solution by the method of lines is then considered, using a finite volume scheme for spatial discretization. An error estimate is derived for the solution of the resulting semidiscrete scheme. Subsequently, sample numerical simulations in two and three dimensions are presented and the experimental convergence measurement is discussed. PubDate: 2017-10-10 DOI: 10.1007/s10543-017-0687-4

Authors:Tomoya Kemmochi Abstract: Abstract In this paper, we develop an energy dissipative numerical scheme for gradient flows of planar curves, such as the curvature flow and the elastic flow. Our study presents a general framework for solving such equations. To discretize the time variable, we use a similar approach to the discrete partial derivative method, which is a structure-preserving method for gradient flows of graphs. For the approximation of curves, we use B-spline curves. Owing to the smoothness of B-spline functions, we can directly address higher order derivatives. Moreover, since B-spline curves require few degrees of freedom, we can reduce the computational cost. In the last part of the paper, we present some numerical examples of the elastic flow, which exhibit topology-changing solutions and more complicated evolution. Videos illustrating our method are available on YouTube. PubDate: 2017-09-27 DOI: 10.1007/s10543-017-0685-6

Authors:Annika Lang; Andreas Petersson; Andreas Thalhammer Abstract: Abstract The (asymptotic) behaviour of the second moment of solutions to stochastic differential equations is treated in mean-square stability analysis. This property is discussed for approximations of infinite-dimensional stochastic differential equations and necessary and sufficient conditions ensuring mean-square stability are given. They are applied to typical discretization schemes such as combinations of spectral Galerkin, finite element, Euler–Maruyama, Milstein, Crank–Nicolson, and forward and backward Euler methods. Furthermore, results on the relation to stability properties of corresponding analytical solutions are provided. Simulations of the stochastic heat equation illustrate the theory. PubDate: 2017-09-14 DOI: 10.1007/s10543-017-0684-7

Authors:Tomasz Hrycak; Sebastian Schmutzhard Abstract: Abstract This paper studies an approximation to the Chebyshev polynomial \(T_n\) computed via a three-term recurrence in floating-point arithmetic. It is shown that close to either endpoint of the interval \([-1, 1]\) , the numerical approximation coincides with the line tangent to \(T_n\) at that endpoint. From this representation new upper and lower error bounds are derived. PubDate: 2017-09-13 DOI: 10.1007/s10543-017-0683-8

Authors:Jing Gao; Arieh Iserles Abstract: Abstract The Filon–Clenshaw–Curtis method (FCC) for the computation of highly oscillatory integrals is known to attain surprisingly high precision. Yet, for large values of frequency \(\omega \) it is not competitive with other versions of the Filon method, which use high derivatives at critical points and exhibit high asymptotic order. In this paper we propose to extend FCC to a new method, FCC \(+\) , which can attain an arbitrarily high asymptotic order while preserving the advantages of FCC. Numerical experiments are provided to illustrate that FCC \(+\) shares the advantages of both familiar Filon methods and FCC, while avoiding their disadvantages. PubDate: 2017-09-12 DOI: 10.1007/s10543-017-0682-9