Search

The matrix-valued Allen-Cahn (MAC) equation was first introduced as a model problem of finding the stationary points of an energy for orthogonal matrix-valued functions and has attracted much attention in recent years. It is well known that the MAC equation satisfies the maximum bound principle (MBP) with respect to either the matrix 2-norm or the Frobenius norm, which plays a key role in understanding the physical meaning and the wellposedness of the model. To preserve this property, we extend the explicit integrating factor Runge-Kutta (IFRK) method to the MAC equation. Moreover, we construct a new three-stage third-order and a new four-stage fourth-order IFRK schemes based on the classical Runge-Kutta schemes. Under a reasonable time-step constraint, we prove the MBP preservation of the IFRK method with respect to the matrix 2-norm, based on which, we further establish their optimal error estimates in the matrix 2-norm. Although numerical results indicate that the IFRK method preserves the MBP with respect to the Frobenius norm, a detailed analysis shows that it is hard to prove this preservation by using the same approach for the case of 2-norm. Several numerical experiments are carried out to test the convergence of the IFRK schemes and to verify the MBP preservation with respect to the matrix 2-norm and the Frobenius norm, respectively. Energy stability is also observed, which clearly indicates the orthogonality of the stationary solution of the MAC equation. In addition, we simulate the coarsening dynamics to verify the motion law of the interface for different initial conditions. [ABSTRACT FROM AUTHOR]

This paper proposes a framework for constructing high-order linearly implicit exponential integrators that conserve a quadratic invariant. This is then applied to the scalar auxiliary variable (SAV) approach. Quadratic invariants are significant objects that are present in various physical equations and also in computationally efficient conservative schemes for general invariants. For instance, the SAV approach converts the invariant into a quadratic form by introducing scalar auxiliary variables, which have been intensively studied in recent years. In this vein, Sato et al. (Appl. Numer. Math. 187, 71-88 2023) proposed high-order linearly implicit schemes that conserve a quadratic invariant. In this study, it is shown that their method can be effectively merged with the Lawson transformation, a technique commonly utilized in the construction of exponential integrators. It is also demonstrated that combining the constructed exponential integrators and the SAV approach yields schemes that are computationally less expensive. Specifically, the main part of the computational cost is the product of several matrix exponentials and vectors, which are parallelizable. Moreover, we conduct some mathematical analyses on the proposed schemes. [ABSTRACT FROM AUTHOR]

In this paper, we introduce a new framework for deriving partitioned implicit-exponential integrators for stiff systems of ordinary differential equations and construct several time integrators of this type. The new approach is suited for solving systems of equations where the forcing term is comprised of several additive nonlinear terms. We analyze the stability, convergence, and efficiency of the new integrators and compare their performance with existing schemes for such systems using several numerical examples. We also propose a novel approach to visualizing the linear stability of the partitioned schemes, which provides a more intuitive way to understand and compare the stability properties of various schemes. Our new integrators are A-stable, second-order methods that require only one call to the linear system solver and one exponential-like matrix function evaluation per time step. [ABSTRACT FROM AUTHOR]

In this paper, we study a new splitting method for the semi-linear system of ordinary differential equation, where the linear part is stiff. Firstly, the stiff part is split into two parts. The first stiff part, that is called the stiff-cutter and expected to be easily inverted, is implicitly treated. The second stiff part and the remaining nonlinear part are explicitly treated. Therefore, such stiff-cut method can be fast implemented and save the CPU time. Theoretically, we rigorously prove that the proposed method is unconditionally stable and convergent, if the stiff-cutter is chosen to be well-matched in the stiff part. As an application, we apply our method to solve a spatial-fractional reaction-diffusion equation and give a way for how to choose a suitable stiff-cutter. Finally, numerical experiments are carried out to illustrate the accuracy and efficiency of the proposed stiff-cut method. [ABSTRACT FROM AUTHOR]

We establish a class of uniformly accurate nested Picard iterative integrator (NPI) Fourier pseudospectral methods for the nonlinear Klein-Gordon-Schrödinger equation (KGS) in the nonrelativistic regime, involving a dimensionless parameter ε ≪ 1 inversely proportional to the speed of light. Actually, the solution propagates waves in time with O(ε2) wavelength when 0 < ε ≪ 1, which brings significant difficulty in designing accurate and efficient numerical schemes. The NPI method is designed by separating the oscillatory part from the non-oscillatory part, and integrating the former exactly. Based on the Picard iteration, the NPI method can be applied to derive arbitrary higher-order methods in time with optimal and uniform accuracy (w.r.t. ε ∈ (0,1]), and the corresponding error estimates are rigorously established. In addition, the practical implementation of the second-order NPI method via Fourier pseupospectral discretization is clearly demonstrated, with extensions to the third-order NPI. Some numerical examples are provided to support our theoretical results and show the accuracy and efficiency of the proposed schemes. [ABSTRACT FROM AUTHOR]

In this study, we introduce a numerical method to reconstruct, from posterior measurements of the solution, the initial condition in the one-dimensional heat conduction problem. The method uses a wavelet-based Galerkin method for the spatial discretization and spectral decomposition of the basis stiffness matrix to get the numerical solution without time discretization as in classical approaches. In fact, according to the proposed method, setting an error bound and properly selecting the eigenvalues of the discretization system, we are able to reconstruct the initial conditions without exceeding the required error. The applicability and computational efficiency of the methods are investigated by solving some numerical examples with toy solutions and the experimental results confirm its accuracy. [ABSTRACT FROM AUTHOR]

In this manuscript, we present a common tensor framework which can be used to generalize one-dimensional numerical tasks to arbitrary dimension d by means of tensor product formulas. This is useful, for example, in the context of multivariate interpolation, multidimensional function approximation using pseudospectral expansions and solution of stiff differential equations on tensor product domains. The key point to obtain an efficient-to-implement BLAS formulation consists in the suitable usage of the μ-mode product (also known as tensor-matrix product or mode-n product) and related operations, such as the Tucker operator. Their MathWorks MATLAB®/GNU Octave implementations are discussed in the paper, and collected in the package KronPACK. We present numerical results on experiments up to dimension six from different fields of numerical analysis, which show the effectiveness of the approach. [ABSTRACT FROM AUTHOR]

Exponential integrators, due to their robust stability properties, have been considered as reliable schemes for numerical solutions of stiff systems. In this paper, we propose generalized exponential time differencing (GETD) schemes for nonlinear fractional differential equations of order α ∈ (0,1). First, we improve the suboptimal performance of the multistep GETD schemes. Using graded mesh, uniform optimal convergence rates under no additional smoothness requirements are obtained. Second, we develop and analyze novel second-order and third-order accurate predictor-corrector type GETD schemes. Using linear stability analysis and numerical illustrations, we demonstrate that the newly introduced schemes have better stability properties than the multistep GETD schemes. Partial fraction decompositions of global Padé approximations for Mittag-Leffler function are used for efficient implementation. Numerical examples involving nonlinear scalar equations and stiff systems are provided to illustrate the theoretical findings. [ABSTRACT FROM AUTHOR]

For the Dirac equation with potentials characterized by a small parameter ε ∈ (0,1], the numerical methods for long-time dynamics have received more and more attention. Recently, two exponential wave integrator Fourier pseudo-spectral (EWIFP) methods for the Dirac equation have been proposed (Feng et al., Appl. Numer. Math. 172, 50–66, 2022) which are uniformly accurate about ε and perform well over the classical methods. However, the EWIFP methods cannot preserve the mass and energy, which are important structural features of the Dirac equation from the perspective of geometric numerical integration. In addition, the EWIFP methods are not time symmetric or only are conditionally stable under specific stability condition which implies CFL condition restrictions on the grid ratio. In this work, we propose a structure-preserving EWIFP (SPEWIFP) method. The proposed method is proved to be time symmetric, stable only under the condition τ ≲ 1 , and preserves the discrete energy and modified mass. Without any CFL condition restrictions on the grid ratio, we carry out a rigourously error analysis and give uniform error bounds of the method at O (h m 0 + ε 1 − β τ 2) up to the time at O(1/εβ) with β ∈ [0,1], mesh size h, time step τ, and an integer m0 determined by the regularity conditions. In general, the Dirac equation with small potentials can be converted to an oscillatory Dirac equation with wavelength at O(εβ) in time which includes the case of simultaneous massless and nonrelativistic regime. It is easy to extend the error bounds and structure-preservation properties to the oscillatory Dirac equation. Numerical experiments support our error bounds and structure-preservation properties. [ABSTRACT FROM AUTHOR]

This paper is concerned with Krylov subspace methods based on the tensor t-product for computing certain quantities associated with generalized third-order tensor functions. We use the tensor t-product and define the tensor global Golub-Kahan bidiagonalization process for approximating tensor functions. Pairs of Gauss and Gauss-Radau quadrature rules are applied to determine the desired quantities with error bounds. An application to the computation of the tensor nuclear norm is presented and illustrates the effectiveness of the proposed methods. [ABSTRACT FROM AUTHOR]

The solution of linear non-autonomous ordinary differential equation systems (also known as the time-ordered exponential) is a computationally challenging problem arising in a variety of applications. In this work, we present and study a new framework for the computation of bilinear forms involving the time-ordered exponential. Such a framework is based on an extension of the non-Hermitian Lanczos algorithm to 4-mode tensors. Detailed results concerning its theoretical properties are presented. Moreover, computational results performed on real-world problems confirm the effectiveness of our approach. [ABSTRACT FROM AUTHOR]

This paper is devoted to the construction and analysis of uniformly accurate (UA) nested Picard iterative integrators (NPI) for highly oscillatory second-order differential equations. The equations involve a dimensionless parameter ε ∈ (0,1], and their solutions are highly oscillatory in time with wavelength at O (ε 2) , which brings severe burdens in numerical computation when ε ≪ 1. In this work, we first propose two NPI schemes for solving a differential equation. The schemes are uniformly first- and second-order accurate for all ε ∈ (0,1]. Moreover, they are super convergent when the time-step size is smaller than ε2. Then, the schemes are generalized to a system of differential equations with the same uniform accuracies. Error bounds are rigorously established and numerical results are reported to confirm the error estimates. [ABSTRACT FROM AUTHOR]

We investigate the problem of approximating the matrix function f(A) by r(A), with f a Markov function, r a rational interpolant of f, and A a symmetric Toeplitz matrix. In a first step, we obtain a new upper bound for the relative interpolation error 1 − r/f on the spectral interval of A. By minimizing this upper bound over all interpolation points, we obtain a new, simple, and sharp a priori bound for the relative interpolation error. We then consider three different approaches of representing and computing the rational interpolant r. Theoretical and numerical evidence is given that any of these methods for a scalar argument allows to achieve high precision, even in the presence of finite precision arithmetic. We finally investigate the problem of efficiently evaluating r(A), where it turns out that the relative error for a matrix argument is only small if we use a partial fraction decomposition for r following Antoulas and Mayo. An important role is played by a new stopping criterion which ensures to automatically find the degree of r leading to a small error, even in presence of finite precision arithmetic. [ABSTRACT FROM AUTHOR]

Let T be a square matrix with a real spectrum, and let f be an analytic function. The problem of the approximate calculation of f(T) is discussed. Applying the Schur triangular decomposition and the reordering, one can assume that T is triangular and its diagonal entries tii are arranged in increasing order. To avoid calculations using the differences tii − tjj with close (including equal) tii and tjj, it is proposed to represent T in a block form and calculate the two main block diagonals using interpolating polynomials. The rest of the f(T) entries can be calculated using the Parlett recurrence algorithm. It is also proposed to perform some scalar operations (such as the building of interpolating polynomials) with an enlarged number of significant decimal digits. [ABSTRACT FROM AUTHOR]

Prior recent work, devoted to the study of polynomial Krylov techniques for the approximation of the action of the matrix exponential etAv, is extended to the case of associated φ-functions (which occur within the class of exponential integrators). In particular, a posteriori error bounds and estimates, based on the notion of the defect (residual) of the Krylov approximation are considered. Computable error bounds and estimates are discussed and analyzed. This includes a new error bound which favorably compares to existing error bounds in specific cases. The accuracy of various error bounds is characterized in relation to corresponding Ritz values of A. Ritz values yield properties of the spectrum of A (specific properties are known a priori, e.g., for Hermitian or skew-Hermitian matrices) in relation to the actual starting vector v and can be computed. This gives theoretical results together with criteria to quantify the achieved accuracy on the fly. For other existing error estimates, the reliability and performance are studied by similar techniques. Effects of finite precision (floating point arithmetic) are also taken into account. [ABSTRACT FROM AUTHOR]

PageRank is a widespread model for analysing the relative relevance of nodes within large graphs arising in several applications. In the current paper, we present a cost-effective Hessenberg-type method built upon the Hessenberg process for the solution of difficult PageRank problems. The new method is very competitive with other popular algorithms in this field, such as Arnoldi-type methods, especially when the damping factor is close to 1 and the dimension of the search subspace is large. The convergence and the complexity of the proposed algorithm are investigated. Numerical experiments are reported to show the efficiency of the new solver for practical PageRank computations. [ABSTRACT FROM AUTHOR]

The concept of Gauss quadrature can be generalized to approximate linear functionals with complex moments. Following the existing literature, this survey will revisit such generalization. It is well known that the (classical) Gauss quadrature for positive definite linear functionals is connected with orthogonal polynomials, and with the (Hermitian) Lanczos algorithm. Analogously, the Gauss quadrature for linear functionals is connected with formal orthogonal polynomials, and with the non-Hermitian Lanczos algorithm with look-ahead strategy; moreover, it is related to the minimal partial realization problem. We will review these connections pointing out the relationships between several results established independently in related contexts. Original proofs of the Mismatch Theorem and of the Matching Moment Property are given by using the properties of formal orthogonal polynomials and the Gauss quadrature for linear functionals. [ABSTRACT FROM AUTHOR]

Two new Runge–Kutta (RK) pairs of orders 6(4) and 7(5) are presented for solving numerically the inhomogeneous linear initial value problems with constant coefficients. These new pairs use only six and eight stages per step respectively. Six stages are needed for conventional Runge–Kutta pairs of orders 5(4) while for such a pair of orders 6(5) we use eight stages. Thus, our proposal is an improvement and it is achieved since the set of order conditions is smaller in the case of interest here. Since traditional simplifications for derivation of Runge–Kutta methods do not apply for this reduced set, we proceed using the differential evolution technique for solving it. We finalize by performing tests over some relevant problems with very pleasant results. [ABSTRACT FROM AUTHOR]

This paper is devoted to the explicit pseudo two-step exponential Runge–Kutta (EPTSERK) methods for the numerical integration of first-order ordinary differential equations. These methods inherit the structure of explicit pseudo two-step Runge–Kutta methods and explicit exponential Runge–Kutta methods. We analyze the order conditions and the global errors of the new methods. The new methods are of order s + 1 with s-stages for some suitable nodes. Numerical experiments are reported to show the convergence and the efficiency of the new methods. [ABSTRACT FROM AUTHOR]

In the last few decades, Runge-Kutta-Nyström (RKN) methods have made significant progress and the study of RKN-type methods for solving highly oscillatory differential equations has received a great deal of attention. In this paper, from the point of view of geometric integration, the oscillation-preserving behaviour of the existing RKN-type methods in the literature are analysed in detail. To this end, it is convenient to introduce the concept of oscillation preservation for RKN-type methods. It turns out that if both the internal stages and the updates of an RKN-type method respect the qualitative and global features of the highly oscillatory solution, then the method is oscillation preserving. Since the internal stages of standard RKN and adapted RKN (ARKN) methods are inimical to the oscillation-preserving structure, neither ARKN methods nor the symplectic and symmetric RKN methods, and standard RKN methods are oscillation preserving. Other concerns relating to oscillation preservation are also considered. In particular, we are concerned with the computational issues for efficiently solving semi-discrete wave equations such as semi-discrete Klein-Gordon (KG) equations and damped sine-Gordon equations. The results of numerical experiments show the importance of the oscillation-preserving property for an RKN-type method and the remarkable superiority of oscillation-preserving integrators when applied to nonlinear multi-frequency highly oscillatory systems. [ABSTRACT FROM AUTHOR]

The main theme of this paper is explicit Gautschi-type integrators for the nonlinear multi-frequency oscillatory second-order initial value problems of the form y ′′ = − A (t , y) y + f (t , y) , y (t 0) = y 0 , y ′ (t 0) = y 0 ′ . This work is important and interesting within the broader framework of the subject. In fact, the Gautschi-type methods for oscillatory problems with a constant matrix A have been investigated by many authors. The key question now is that the classical variation-of-constants approach is not applicable to the oscillatory nonlinear problems with a variable coefficient matrix A(t,y). We consider successive approximations or locally equivalent systems for the problems, and derive efficient explicit Gautschi-type integrators. The error analysis is presented for the local approximation accordingly. Accompanying numerical results demonstrate the remarkable efficiency of the new Gautschi-type integrators in comparison with some existing numerical methods in the scientific literature. [ABSTRACT FROM AUTHOR]

This paper analyzes global error bounds of one-stage explicit extended Runge–Kutta–Nyström integrators for semilinear wave equations. The analysis is presented by using spatial semidiscretizations with periodic boundary conditions. Optimal second-order convergence is proved without requiring Lipschitz continuous and higher regularity of the exact solution. Moreover, the error analysis is not restricted to the spectral semidiscretization in space. [ABSTRACT FROM AUTHOR]

Multi-frequency, highly oscillatory Hamiltonian problems derive from the mathematical modelling of many real-life applications. We here propose a variant of Hamiltonian Boundary Value Methods (HBVMs), which is able to efficiently deal with the numerical solution of such problems. We present algorithms to select the parameters of the methods that allow one to obtain numerical approximations with spectral accuracy. We also propose an efficient implementation of the methods when using a Newton-type iteration to solve the implicit equations associated with this class of formulas. [ABSTRACT FROM AUTHOR]

In last few years, many ERKN methods have been investigated for solving multi-frequency multidimensional second-order ordinary differential equations, and the numerical efficiency has been checked strongly in scientific computation. But in the constructions of (especially high-order) new ERKN methods, lots of time and effort are costed in presenting the practical order conditions firstly and then in adding some reasonable assumptions to get the coefficient functions finally. In this paper, a feasible and effective technique is given which makes the construction of ERKN methods finished in a few seconds or a few minutes, even for high-order integrators. Moreover, this technique does not need any more information and knowledge except the classical RKN method. And this paper also gives the theoretical explanation to guarantee that the ERKN method obtained from this technique has the same order and the same properties as the underlying RKN method. [ABSTRACT FROM AUTHOR]

We present an extended Krylov subspace analogue of the two-sided Lanczos method, i.e., a method which, given a nonsingular matrix A and vectors b, c with $\left \langle {{\mathbf {b}},{\mathbf {c}}}\right \rangle \neq 0$ , constructs bi-orthonormal bases of the extended Krylov subspaces ${\mathcal {E}}_{m}(A,{\mathbf {b}})$ and ${\mathcal {E}}_{m}(A^{T}\!,{\mathbf {c}})$ via short recurrences. We investigate the connection of the proposed method to rational moment matching for bilinear forms c f( A) b, similar to known results connecting the two-sided Lanczos method to moment matching. Numerical experiments demonstrate the quality of the resulting approximations and the numerical behavior of the new extended Krylov subspace method. [ABSTRACT FROM AUTHOR]

We assess the applicability and efficiency of time-adaptive high-order splitting methods applied for the numerical solution of (systems of) nonlinear parabolic problems under periodic boundary conditions. We discuss in particular several applications generating intricate patterns and displaying nonsmooth solution dynamics. First, we give a general error analysis for splitting methods for parabolic problems under periodic boundary conditions and derive the necessary smoothness requirements on the exact solution in particular for the Gray-Scott equation and the Van der Pol equation. Numerical examples demonstrate the convergence of the methods and serve to compare the efficiency of different time-adaptive splitting schemes and of splitting into either two or three operators, based on appropriately constructed a posteriori local error estimators. [ABSTRACT FROM AUTHOR]

A fourth-order implicit-explicit time-discretization scheme based on the exponential time differencing approach with a fourth-order compact scheme in space is proposed for space fractional nonlinear Schrödinger equations. The stability and convergence of the compact scheme are discussed. It is shown that the compact scheme is fourth-order convergent in space and in time. Numerical experiments are performed on single and coupled systems of two and four fractional nonlinear Schrödinger equations. The results demonstrate accuracy, efficiency, and reliability of the scheme. A linearly implicit conservative method with the fourth-order compact scheme in space is also considered and used on the system of space fractional nonlinear Schrödinger equations. [ABSTRACT FROM AUTHOR]

In this paper, we present a high-order accurate method for two-dimensional semilinear parabolic equations. The method is based on a Galerkin-Chebyshev spectral method for discretizing spatial derivatives and a block boundary value methods of fourth-order for temporal discretization. Our formulation has high-order accurate in both space and time. Optimal a priori error bound is derived in the weighted $L^{2}_{\omega }$-norm for the semidiscrete formulation. Extensive numerical results are presented to demonstrate the convergence properties of the method. [ABSTRACT FROM AUTHOR]

New algorithms are developed for estimating the condition number of f( A) b, where A is a matrix and b is a vector. The condition number estimation algorithms for f( A) already available in the literature require the explicit computation of matrix functions and their Fr´echet derivatives and are therefore unsuitable for the large, sparse A typically encountered in f( A) b problems. The algorithms we propose here use only matrix-vector multiplications. They are based on a modified version of the power iteration for estimating the norm of the Fr´echet derivative of a matrix function, and work in conjunction with any existing algorithm for computing f( A) b. The number of matrix-vector multiplications required to estimate the condition number is proportional to the square of the number of matrix-vector multiplications required by the underlying f( A) b algorithm. We develop a specific version of our algorithm for estimating the condition number of e b, based on the algorithm of Al-Mohy and Higham (SIAM J. Matrix Anal. Appl. 30(4), 1639-1657, ). Numerical experiments demonstrate that our condition estimates are reliable and of reasonable cost. [ABSTRACT FROM AUTHOR]

We consider the numerical integration of linear-quadratic optimal control problems. This problem requires the solution of a boundary value problem: a non-autonomous matrix Riccati differential equation (RDE) with final conditions coupled with the state vector equation with initial conditions. The RDE has positive definite matrix solution and to numerically preserve this qualitative property we propose first to integrate this equation backward in time with a sufficiently accurate scheme. Then, this problem turns into an initial value problem, and we analyse splitting and Magnus integrators for the forward time integration which preserve the positive definite matrix solutions for the RDE. Duplicating the time as two new coordinates and using appropriate splitting methods, high order methods preserving the desired property can be obtained. The schemes make sequential computations and do not require the storrage of intermediate results, so the storage requirements are minimal. The proposed methods are also adapted for solving linear-quadratic N-player differential games. The performance of the splitting methods can be considerably improved if the system is a perturbation of an exactly solvable problem and the system is properly split. Some numerical examples illustrate the performance of the proposed methods. [ABSTRACT FROM AUTHOR]

Krylov subspace methods for approximating a matrix function f( A) times a vector v are analyzed in this paper. For the Arnoldi approximation to e v, two reliable a posteriori error estimates are derived from the new bounds and generalized error expansion we establish. One of them is similar to the residual norm of an approximate solution of the linear system, and the other one is determined critically by the first term of the error expansion of the Arnoldi approximation to e v due to Saad. We prove that each of the two estimates is reliable to measure the true error norm, and the second one theoretically justifies an empirical claim by Saad. In the paper, by introducing certain functions ϕ( z) defined recursively by the given function f( z) for certain nodes, we obtain the error expansion of the Krylov-like approximation for f( z) sufficiently smooth, which generalizes Saad's result on the Arnoldi approximation to e v. Similarly, it is shown that the first term of the generalized error expansion can be used as a reliable a posteriori estimate for the Krylov-like approximation to some other matrix functions times v. Numerical examples are reported to demonstrate the effectiveness of the a posteriori error estimates for the Krylov-like approximations to e v, cos( A) v and sin( A) v. [ABSTRACT FROM AUTHOR]

In this paper, we are concerned with the error analysis for the two-step extended Runge-Kutta-Nyström-type (TSERKN) methods [Comput. Phys. Comm. 182 (2011) 2486-2507] for multi-frequency and multidimensional oscillatory systems y″( t) + My( t) = f( t, y( t)), where high-frequency oscillations in the solutions are generated by the linear part My( t). TSERKN methods extend the two-step hybrid methods [IMA J. Numer. Anal. 23 (2003) 197-220] by reforming both the internal stages and the updates so that they are adapted to the oscillatory properties of the exact solutions. However, the global error analysis for the TSERKN methods has not been investigated. In this paper we construct a new three-stage explicit TSERKN method of order four and present the global error bound for the new method, which is proved to be independent of ∥ M∥ under suitable assumptions. This property of our new method is very important for solving highly oscillatory systems (1), where ∥ M∥ may be arbitrarily large. We also analyze the stability and phase properties for the new method. Numerical experiments are included and the numerical results show that the new method is very competitive and promising compared with the well-known high quality methods proposed in the scientific literature. [ABSTRACT FROM AUTHOR]

The numerical integration of Hamiltonian systems with multi-frequency and multidimensional oscillatory solutions is encountered frequently in many fields of the applied sciences. In this paper, we firstly summarize the extended Runge-Kutta-Nyström (ERKN) methods proposed by Wu et al. (Comput. Phys. Comm. 181:1873-1887, (2010)) for multi-frequency and multidimensional oscillatory systems and restate the order conditions and symplecticity conditions for the explicit ERKN methods. Secondly, we devote to exploring the explicit symplectic multi-frequency and multidimensional ERKN methods of order five based on the symplecticity conditions and order conditions. A five-stage explicit symplectic multi-frequency and multidimensional ERKN method of order five with some small residuals is proposed and its stability and phase properties are analyzed. It is shown that the new method is dispersive of order six. Numerical experiments are carried out and the numerical results demonstrate that the new method is much more efficient than the methods appeared in the scientific literature. [ABSTRACT FROM AUTHOR]

We are interested in computing the Fermi-Dirac matrix function in which the matrix argument is the Hamiltonian matrix arising from density functional theory (DFT) applications. More precisely, we are really interested in the diagonal of this matrix function. We discuss rational approximation methods to the problem, specifically the rational Chebyshev approximation and the continued fraction representation. These schemes are further decomposed into their partial fraction expansions, leading ultimately to computing the diagonal of the inverse of a shifted matrix over a series of shifts. We describe Lanczos and sparse direct methods to address these systems. Each approach has advantages and disadvantages that are illustrated with experiments. [ABSTRACT FROM AUTHOR]

Abstract  Given a nonsingular complex matrix and complex vectors v and w of length N, one may wish to estimate the quadratic form w * A  − 1 v, where w * denotes the conjugate transpose of w. This problem appears in many applications, and Gene Golub was the key figure in its investigations for decades. He focused mainly on the case A Hermitian positive definite (HPD) and emphasized the relationship of the algebraically formulated problems with classical topics in analysis - moments, orthogonal polynomials and quadrature. The essence of his view can be found in his contribution Matrix Computations and the Theory of Moments, given at the International Congress of Mathematicians in Zürich in 1994. As in many other areas, Gene Golub has inspired a long list of coauthors for work on the problem, and our contribution can also be seen as a consequence of his lasting inspiration. In this paper we will consider a general mathematical concept of matching moments model reduction, which as well as its use in many other applications, is the basis for the development of various approaches for estimation of the quadratic form above. The idea of model reduction via matching moments is well known and widely used in approximation of dynamical systems, but it goes back to Stieltjes, with some preceding work done by Chebyshev and Heine. The algebraic moment matching problem can for A HPD be formulated as a variant of the Stieltjes moment problem, and can be solved using Gauss-Christoffel quadrature. Using the operator moment problem suggested by Vorobyev, we will generalize model reduction based on matching moments to the non-Hermitian case in a straightforward way. Unlike in the model reduction literature, the presented proofs follow directly from the construction of the Vorobyev moment problem. [ABSTRACT FROM AUTHOR]