Search

This paper considers computing partial eigenpairs of differential eigenvalue problems (DEPs) such that eigenvalues are in a certain region on the complex plane. Recently, based on a "solve-then-discretize" paradigm, an operator analogue of the FEAST method has been proposed for DEPs without discretization of the coefficient operators. Compared to conventional "discretize-then-solve" approaches that discretize the operators and solve the resulting matrix problem, the operator analogue of FEAST exhibits much higher accuracy; however, it involves solving a large number of ordinary differential equations (ODEs). In this paper, to reduce the computational costs, we propose operation analogues of Sakurai-Sugiura-type complex moment-based eigensolvers for DEPs using higher-order complex moments and analyze the error bound of the proposed methods. We show that the number of ODEs to be solved can be reduced by a factor of the degree of complex moments without degrading accuracy, which is verified by numerical results. Numerical results demonstrate that the proposed methods are over five times faster compared with the operator analogue of FEAST for several DEPs while maintaining almost the same high accuracy. This study is expected to promote the "solve-then-discretize" paradigm for solving DEPs and contribute to faster and more accurate solutions in real-world applications., Comment: 26 pages, 9 figures

In many applications, it is useful to use piecewise polynomials that satisfy certain regularity conditions at the joint points. Cubic spline functions emerge as good candidates having C2 regularity. On the other hand, if the data points present discontinuities, the classical spline approximations produce Gibbs oscillations. In a recent paper, we have introduced a new nonlinear spline approximation avoiding the presence of these oscillations. Unfortunately, this new reconstruction loses the C2 regularity. This paper introduces a new nonlinear spline that preserves the regularity at all the joint points except at the end points of an interval containing a discontinuity, and that avoids the Gibbs oscillations. Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. This work was funded by the Programa de Apoyo a la investigación de la fundación Séneca-Agencia de Ciencia y Tecnología de la Región de Murcia 20928/PI/18, by the national research project MMTM2015-64382-P and PID2019-108336GB-I00 (MINECO/FEDER), by grant MTM2017-83942 funded by Spanish MINECO and by grant PID2020-117211GB-I00 funded by MCIN/AEI/10.13039/501100011033.

The approximation of trace(f(Ω)), where f is a function of a symmetric matrix Ω, can be challenging when Ω is exceedingly large. In such a case even the partial Lanczos decomposition of Ω is computationally demanding and the stochastic method investigated by Bai et al. (J. Comput. Appl. Math. 74:71–89, 1996) is preferred. Moreover, in the last years, a partial global Lanczos method has been shown to reduce CPU time with respect to partial Lanczos decomposition. In this paper we review these techniques, treating them under the unifying theory of measure theory and Gaussian integration. This allows generalizing the stochastic approach, proposing a block version that collects a set of random vectors in a rectangular matrix, in a similar fashion to the partial global Lanczos method. We show that the results of this technique converge quickly to the same approximation provided by Bai et al. (J. Comput. Appl. Math. 74:71–89, 1996), while the block approach can leverage the same computational advantages as the partial global Lanczos. Numerical results for the computation of the Von Neumann entropy of complex networks prove the robustness and efficiency of the proposed block stochastic method.

In this paper, we analyze the drift-implicit (or backward) Euler numerical scheme for a class of stochastic differential equations with unbounded drift driven by an arbitrary $λ$-Hölder continuous process, $λ\in(0,1)$. We prove that, under some mild moment assumptions on the Hölder constant of the noise, the $L^r(Ω;L^\infty([0,T]))$-rate of convergence is equal to $λ$. To exemplify, we consider numerical schemes for the generalized Cox--Ingersoll-Ross and Tsallis--Stariolo--Borland models. The results are illustrated by simulations., 24 pages, 3 figures

The analysis of linear ill-posed problems often is carried out in function spaces using tools from functional analysis. However, the numerical solution of these problems typically is computed by first discretizing the problem and then applying tools from finite-dimensional linear algebra. The present paper explores the feasibility of applying the Chebfun package to solve ill-posed problems with a regularize-first approach numerically. This allows a user to work with functions instead of vectors and with integral operators instead of matrices. The solution process therefore is much closer to the analysis of ill-posed problems than standard linear algebra-based solution methods. Furthermore, the difficult process of explicitly choosing a suitable discretization is not required.

This paper is concerned with the inverse problem of constructing a symmetric nonnegative matrix from realizable spectrum. We reformulate the inverse problem as an underdetermined nonlinear matrix equation over a Riemannian product manifold. To solve it, we develop a Riemannian underdetermined inexact Newton dogleg method for solving a general underdetermined nonlinear equation defined between Riemannian manifolds and Euclidean spaces. The global and quadratic convergence of the proposed method is established under some mild assumptions. Then we solve the inverse problem by applying the proposed method to its equivalent nonlinear matrix equation and a preconditioner for the perturbed normal Riemannian Newton equation is also constructed. Numerical tests show the efficiency of the proposed method for solving the inverse problem., 32 pages, 6 figures

In this paper, we focus on developing a high efficient algorithm for solving d-dimension time-fractional diffusion equation (TFDE). For TFDE, the initial function or source term is usually not smooth, which can lead to the low regularity of exact solution. And such low regularity have a marked impact on the convergence rate of numerical method. In order to improve the convergence rate of the algorithm, we introduce the space-time sparse grid (STSG) method to solve TFDE. In our study, we employ the sine basis for spatial discretization, and all the sine coefficients can be divided into several levels. The sine coefficients with different levels are discretized by temporal basis with different scales, which can lead to the STSG method. Under certain conditions, the function approximation on standard STSG can achieve the accuracy order O(2-Jd) with O(2JJ) degrees of freedom (DOF) for d=1 and O(2Jd) DOF for d>1, where J denotes the maximal level of sine coefficients. However, the standard STSG is not suitable to simulate the singularity of TFDE at the initial time. To overcome this, we integrate the full grid into the STSG, and obtain the modified STSG. Then, the modified STSG is used to construct the fully discrete scheme for solving TFDE. The great advantage of STSG method can be shown in the comparative numerical experiment.Mathematics Subject Classification (2010) 35R11 · 65M70 · 65T40 · 68Q25

Consider the problem of solving systems of linear algebraic equations $$\varvec{A}\varvec{x}=\varvec{b}$$ A x = b with a real symmetric positive definite matrix $$\varvec{A}$$ A using the conjugate gradient (CG) method. To stop the algorithm at the appropriate moment, it is important to monitor the quality of the approximate solution. One of the most relevant quantities for measuring the quality of the approximate solution is the $$\varvec{A}$$ A -norm of the error. This quantity cannot be easily computed; however, it can be estimated. In this paper we discuss and analyze the behavior of the Gauss-Radau upper bound on the $$\varvec{A}$$ A -norm of the error, based on viewing CG as a procedure for approximating a certain Riemann-Stieltjes integral. This upper bound depends on a prescribed underestimate $$\varvec{\mu }$$ μ to the smallest eigenvalue of $$\varvec{A}$$ A . We concentrate on explaining a phenomenon observed during computations showing that, in later CG iterations, the upper bound loses its accuracy, and is almost independent of $$\varvec{\mu }$$ μ . We construct a model problem that is used to demonstrate and study the behavior of the upper bound in dependence of $$\varvec{\mu }$$ μ , and developed formulas that are helpful in understanding this behavior. We show that the above-mentioned phenomenon is closely related to the convergence of the smallest Ritz value to the smallest eigenvalue of $$\varvec{A}$$ A . It occurs when the smallest Ritz value is a better approximation to the smallest eigenvalue than the prescribed underestimate $$\varvec{\mu }$$ μ . We also suggest an adaptive strategy for improving the accuracy of the upper bounds in the previous iterations.

Algebraically stabilized finite element discretizations of scalar steady-state convection–diffusion–reaction equations often provide accurate approximate solutions satisfying the discrete maximum principle (DMP). However, it was observed that a deterioration of the accuracy and convergence rates may occur for some problems if meshes without local symmetries are used. The paper investigates these phenomena both numerically and analytically and the findings are used to design a new algebraic stabilization called Symmetrized Monotone Upwind-type Algebraically Stabilized (SMUAS) method. It is proved that the SMUAS method is linearity preserving and satisfies the DMP on arbitrary simplicial meshes. Moreover, numerical results indicate that the SMUAS method leads to optimal convergence rates on general simplicial meshes.

In this paper, we consider the Gauss-Kronrod quadrature formulas for a modified Chebyshev weight. Efficient estimates of the error of these Gauss–Kronrod formulae for analytic functions are obtained, using techniques of contour integration that were introduced by Gautschi and Varga (cf. Gautschi and Varga SIAM J. Numer. Anal. 20, 1170–1186 1983). Some illustrative numerical examples which show both the accuracy of the Gauss–Kronrod formulas and the sharpness of our estimations are displayed. Though for the sake of brevity we restrict ourselves to the first kind Chebyshev weight, a similar analysis may be carried out for the other three Chebyshev type weights; part of the corresponding computations are included in a final appendix.

Recently, the efficient numerical solution of Hamiltonian problems has been tackled by defining the class of energy-conserving Runge-Kutta methods named Hamiltonian Boundary Value Methods (HBVMs). Their derivation relies on the expansion of the vector field along the Legendre orthonormal basis. Interestingly, this approach can be extended to cope with other orthonormal bases and, in particular, we here consider the case of the Chebyshev polynomial basis. The corresponding Runge-Kutta methods were previously obtained by Costabile and Napoli [33]. In this paper, the use of a different framework allows us to carry out a novel analysis of the methods also when they are used as spectral formulae in time, along with some generalizations of the methods., Comment: 25 pages, 2 figures, 2 tables

Spectral conjugate gradient (SCG) methods are combinations of spectral gradient method and conjugate gradient (CG) methods, which have been well studied in Euclidean space. In this paper, we aim to extend this class of methods to solve optimization problems on Riemannian manifolds. Firstly, we present a Riemannian version of the spectral parameter, which guarantees that the search direction always satisfies the sufficient descent property without the help of any line search strategy. Secondly, we introduce a generic algorithmic framework for the Riemannian SCG methods, in which the selection of the CG parameter is very flexible. Under the Riemannian Wolfe conditions, the global convergence of the proposed algorithmic framework is established whenever the absolute value of the CG parameter is no more than the Riemannian Fletcher-Reeves CG parameter. Finally, some preliminary numerical results are reported and compared with several classical Riemannian CG methods, which show that our new methods are efficient.MSC Classification: 65K05 , 90C30

In a recent paper, a conformable fractional Newton-type method was proposed for solving nonlinear equations. This method involves a lower computational cost compared to others fractional iterative methods. Indeed, the theoretical order of convergence is held in practice, and it presents a better numerical behaviour than fractional Newton-type methods formerly proposed, even compared to classical Newton-Raphson method. In this work, we design a generalization of this method for solving nonlinear systems by using a new conformable fractional Jacobian matrix, and a suitable conformable Taylor power series; and it is compared with classical Newton’s scheme. The necessary concepts and results are stated in order to design this method. Convergence analysis is made and a quadratic order of convergence is obtained, as in classical Newton's method. Numerical tests are made, and the Approximated Computational Order ofo Convergence (ACOC) supports the theory. Also, the proposed scheme shows good stability properties observed by means of convergence planes.

In this paper, the solution to a bivariate Appell interpolation problem proposed in a previous work is given. Bounds of the truncation error are considered. Ten new interpolants for real, regular, bivariate functions are constructed. Numerical examples and comparisons with bivariate Bernstein polynomials are considered.

Two-level Schwarz domain decomposition methods are very powerful techniques for the efficient numerical solution of partial differential equations (PDEs). A two-level domain decomposition method requires two main components: a one-level preconditioner (or its corresponding smoothing iterative method), which is based on domain decomposition techniques, and a coarse correction step, which relies on a coarse space. The coarse space must properly represent the error components that the chosen one-level method is not capable to deal with. In the literature, most of the works introduced efficient coarse spaces obtained as the span of functions defined on the entire space domain of the considered PDE. Therefore, the corresponding two-level preconditioners and iterative methods are defined in volume. In this paper, we use the excellent smoothing properties of Schwarz domain decomposition methods to define, for general elliptic problems, a new class of substructured two-level methods, for which both Schwarz smoothers and coarse correction steps are defined on the interfaces (except for the application of the smoother that requires volumetric subdomain solves). This approach has several advantages. On the one hand, the required computational effort is cheaper than the one required by classical volumetric two-level methods. On the other hand, our approach does not require, like classical multi-grid methods, the explicit construction of coarse spaces, and it permits a multilevel extension, which is desirable when the high dimension of the problem or the scarce quality of the coarse space prevents the efficient numerical solution. Numerical experiments demonstrate the effectiveness of the proposed new numerical framework.

Null-space methods have long been used to solve large sparse n × n symmetric saddle point systems of equations in which the (2, 2) block is zero. This paper focuses on the case where the (1, 1) block is ill conditioned or rank deficient and the k × k (2, 2) block is non zero and small (k ≪ n). Additionally, the (2, 1) block may be rank deficient. Such systems arise in a range of practical applications. A novel null-space approach is proposed that transforms the system matrix into a nicer symmetric saddle point matrix of order n that has a non zero (2, 2) block of order at most 2k and, importantly, the (1, 1) block is symmetric positive definite. Success of any null-space approach depends on constructing a suitable null-space basis. We propose methods for wide matrices having far fewer rows than columns with the aim of balancing stability of the transformed saddle point matrix with preserving sparsity in the (1, 1) block. Linear least squares problems that contain a small number of dense rows are an important motivation and are used to illustrate our ideas and to explore their potential for solving large-scale systems.

In this paper, we present a novel strategy to systematically construct linearly implicit energy-preserving schemes with arbitrary order of accuracy for Hamiltonian PDEs. Such novel strategy is based on the newly developed exponential scalar variable (ESAV) approach that can remove the bounded-from-blew restriction of nonlinear terms in the Hamiltonian functional and provides a totally explicit discretization of the auxiliary variable without computing extra inner products, which make it more effective and applicable than the traditional scalar auxiliary variable (SAV) approach. To achieve arbitrary high-order accuracy and energy preservation, we utilize the symplectic Runge-Kutta method for both solution variables and the auxiliary variable, where the values of internal stages in nonlinear terms are explicitly derived via an extrapolation from numerical solutions already obtained in the preceding calculation. A prediction-correction strategy is proposed to further improve the accuracy. Fourier pseudo-spectral method is then employed to obtain fully discrete schemes. Compared with the SAV schemes, the solution variables and the auxiliary variable in these ESAV schemes are now decoupled. Moreover, when the linear terms are of constant coefficients, the solution variables can be explicitly solved by using the fast Fourier transform. Numerical experiments are carried out for three Hamiltonian PDEs to demonstrate the efficiency and conservation of the ESAV schemes.

In this paper, we consider a regularized least squares problem subject to convex constraints. Our algorithm is based on the superiorization technique, equipped with a new step size rule which uses subgradient projections. The superiorization method is a two-step method where one step reduces the value of the penalty term and the other step reduces the residual of the underlying linear system (using an algorithmic operator T). For the new step size rule, we present a convergence analysis for the case when T belongs to a large subclass of strictly quasi-nonexpansive operators. To examine our algorithm numerically, we consider box constraints and use the total variation (TV) functional as a regularization term. The specific test cases are chosen from computed tomography using both noisy and noiseless data. We compare our algorithm with previously used parameters in superiorization. The T operator is based on sequential block iteration (for which our convergence analysis is valid), but we also use the conjugate gradient method (without theoretical support). Finally, we compare with the well-known "fast iterative shrinkage-thresholding algorithm" (FISTA). The numerical results demonstrate that our new step size rule improves previous step size rules for the superiorization methodology and is competitive with, and in several instances behaves better than, the other methods.

In the present paper, we consider the semilocal convergence issue of two-step Newton method for solving nonlinear operator equation in Banach spaces. Under the assumption that the first derivative of the operator satisfies a generalized Lipschitz condition, a new semilocal convergence analysis for the two-step Newton method is presented. The Q-cubic convergence is obtained by an additional condition. This analysis also allows us to obtain three important spacial cases about the convergence results based on the premises of Kantorovich, Smale and Nesterov-Nemirovskii types. As applications of our convergence results, a nonsymmetric algebraic Riccati equation arising from transport theory and a two-dimensional nonlinear convection-diffusion equation are provided.

When exploring the literature, it can be observed that the operator obtained when applying Newton-like root finding algorithms to the quadratic polynomials z2 − c has the same form regardless of which algorithm has been used. In this paper, we justify why this expression is obtained. This is done by studying the symmetries of the operators obtained after applying Newton-like algorithms to a family of degree d polynomials p(z) = zd − c. Moreover, we provide an iterative procedure to obtain the expression of new Newton-like algorithms. We also carry out a dynamical study of the given generic operator and provide general conclusions of this type of methods.

This paper concerns test matrices for numerical linear algebra using an error-free transformation of floating-point arithmetic. For specified eigenvalues given by a user, we propose methods of generating a matrix whose eigenvalues are exactly known based on, for example, Schur or Jordan normal form and a block diagonal form. It is also possible to produce a real matrix with specified complex eigenvalues. Such test matrices with exactly known eigenvalues are useful for numerical algorithms in checking the accuracy of computed results. In particular, exact errors of eigenvalues can be monitored. To generate test matrices, we first propose an error-free transformation for the product of three matrices YSX. We approximate S by ${S^{\prime }}$ S ′ to compute ${YS^{\prime }X}$ Y S ′ X without a rounding error. Next, the error-free transformation is applied to the generation of test matrices with exactly known eigenvalues. Note that the exactly known eigenvalues of the constructed matrix may differ from the anticipated given eigenvalues. Finally, numerical examples are introduced in checking the accuracy of numerical computations for symmetric and unsymmetric eigenvalue problems.