Search

In this paper, we survey selected software packages for the numerical solution of boundary value ODEs (BVODEs), time-dependent PDEs in one spatial dimension (1DPDEs), and initial value ODEs (IVODEs). A unifying theme of this paper is our focus on software packages for these problem classes that compute error-controlled, continuous numerical solutions. A continuous numerical solution can be accessed by the user at any point in the domain. We focus on error-control software; this means that the software adapts the computation until it obtains a continuous approximate solution with a corresponding error estimate that satisfies the user tolerance. The second section of the paper will provide an overview of recent work on the development of COLNEWSC, an updated version of the widely used collocation BVODE solver, COLNEW, that returns an error-controlled continuous approximate solution based on the use of a superconvergent interpolant to the underlying collocation solution. The third section of the paper gives a brief review of recent work on the development of a new 1DPDE solver, BACOLIKR, that provides time- and space-dependent event detection for an error-controlled continuous numerical solution. In the fourth section of the paper, we briefly review the state of the art in IVODE software for the computation of error-controlled continuous numerical solutions. [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]

Runge–Kutta embedded pairs of high algebraic order are frequently utilized when strict tolerances are required. When creating such pairings of orders nine and eight for use in double precision arithmetic, numerous conditions are often satisfied. First and foremost, we strive to keep the coefficients' magnitudes small to prevent accuracy loss. We may, however, allow greater coefficients when working with quadruple precision. Then, we may build pairs of orders 9 and 8 with significantly smaller truncation errors. In this paper, a novel pair is generated that, as predicted, outperforms state-of-the-art pairs of the same orders in a collection of important problems. [ABSTRACT FROM AUTHOR]

This paper develops an improved tri-coloured rooted-tree theory for the order conditions for ERKN methods solving general multi-frequency and multidimensional second-order oscillatory systems. The bottleneck of the original tricoloured rooted-tree theory is the existence of numerous redundant trees. In light of the fact that the sum of the products of the symmetries and the elementary differentials is meaningful, this paper naturally introduces the so-called extended elementary differential mappings. Then, the new improved tri-coloured rooted tree theory is established based on a subset of the original tri-coloured rooted-tree set. This new theory makes all redundant trees disappear, and thus, the order conditions of ERKN methods for general multi-frequency and multidimensional second-order oscillatory systems are reduced greatly. Furthermore, with this new theory, we present some new ERKN methods of order up to four. Numerical experiments are implemented and the results show that ERKN methods can be competitive with other existing methods in the scientific literature, especially when comparatively large stepsizes are used. [ABSTRACT FROM AUTHOR]

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 (Numer. Algo., 27, 119–130 2021). 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. [ABSTRACT FROM AUTHOR]

In this paper, the system of nonlinear inequalities is considered. The problem is approximated by the parameterized smooth equations which is formed by Chen-Harker-Kanzow-Smale smoothing function. A smoothing homotopy method is proposed to solve the system of smooth equations. For almost all initial points in Euclidean space, the existence and convergence of the homotopy path is proved under certain weaker condition without assumption of the monotonicity. The numerical results illustrate that the method is effective. [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]

The connection between trees and differential equations was pointed out in the classic paper by Cayley (Phil. Mag. 13, 172–176 1857). Trees were also used in the work of Merson (1957), on the order of Runge–Kutta methods. The paper by Hairer and Wanner (Computing 13, 1–15 1974), where the term B-series was introduced, followed papers by the present author (Butcher J. Austral. Math. Soc. 3, 185–201 1963, Math. Comput. 26, 79–106 1972). The present paper will survey the use of trees in the formulation of B-series and illustrate the results by constructing and analysing some examples of general linear methods. [ABSTRACT FROM AUTHOR]

In this paper, we are concerned with energy-conserving methods for Poisson problems, which are effectively solved by defining a suitable generalization of HBVMs, a class of energy-conserving methods for Hamiltonian problems. The actual implementation of the methods is fully discussed, with a particular emphasis on the conservation of Casimirs. Some numerical tests are reported, in order to assess the theoretical findings. [ABSTRACT FROM AUTHOR]

This paper is devoted to investigate the modified extended second derivative backward differentiation formulae from second derivative general linear methods point of view. This makes it possible to open some maneuver rooms in developing the methods with superior features by perturbing the abscissa vector of the methods. The proposed methods are constructed to have better accuracy and stability properties in comparison with the original ones. These improvements are verified by giving some numerical experiments. [ABSTRACT FROM AUTHOR]

In this paper, we describe an approach to estimate the global error for explicit second derivative general linear methods based on the approach which has been already used for global error estimation of explicit general linear methods. In this approach, to estimate the global error, we use the numerical solutions of pairs of second derivative general linear methods with the same order and stage order that are constructed such that their global error functions are proportional. Numerical experiments demonstrate the excellent agreement of the global error estimation with the exact one in both constant and variable stepsize environments. [ABSTRACT FROM AUTHOR]

In this paper, we exploit the Strang splitting technique for solving the multidimensional Allen-Cahn equations with anisotropic spatial fractional Riesz derivatives. Fully discrete numerical methods are proposed using exponential Strang's splitting schemes for the time integration with finite difference discretization in space. It is proved that the proposed methods can preserve the discrete maximum principle unconditionally. Furthermore, the fully discrete methods are theoretically confirmed to be convergent with second-order accuracy in both of time and space. In practical implementation, the proposed algorithms require to compute the matrix exponential associated with only one-dimensional discretized matrices that possess Toeplitz structure. Meanwhile, a fast algorithm is further developed for evaluating the product of the Toeplitz matrix exponential with a vector. Numerical examples are presented to verify the theoretical analysis and demonstrate the efficiency of the proposed methods. [ABSTRACT FROM AUTHOR]

Because the expense of estimating the optimal value of the relaxation parameter in the successive over-relaxation (SOR) method is usually prohibitive, the parameter is often adaptively controlled. In this paper, new adaptive SOR methods are presented that are applicable to a variety of symmetric positive definite linear systems and do not require additional matrix-vector products when updating the parameter. To this end, we regard the SOR method as an algorithm for minimising a certain objective function, which yields an interpretation of the relaxation parameter as the step size following a certain change of variables. This interpretation enables us to adaptively control the step size based on some line search techniques, such as the Wolfe conditions. Numerical examples demonstrate the favourable behaviour of the proposed methods. [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 concerns the discrete time waveform relaxation (DWR) methods for ordinary differential equations (ODEs). We present a general algorithm of constructing the DWR methods with any order of convergence, which applies any numerical methods of ODEs to the perturbed equations of iterative schemes of continuous time waveform relaxation methods. It is demonstrated that the DWR method presented in this paper has the same convergent order as the numerical method used to discretize perturbed equations. Two classes of interpolation polynomials are given to generate perturbed equations. Finally, numerical experiments are presented in order to check against results obtained. [ABSTRACT FROM AUTHOR]

For their several attractive features from the viewpoint of the numerical computations, linear barycentric rational interpolants have been recently used to construct various numerical methods for solving different classes of equations. In this paper, we introduce a family of linear multistep second derivative methods together with a starting procedure based on barycentric rational interpolants. The order of convergence and linear stability properties of the proposed methods have been investigated. To validate the theoretical results and efficiency of the methods, some numerical experiments have been provided. [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]

General Linear Methods (GLMs) were introduced as the natural generalizations of the classical Runge-Kutta and linear multistep methods. An extension of GLMs, so-called SGLMs (GLM with second derivative), was introduced to the case in which second derivatives, as well as first derivatives, can be calculated. In this paper, we introduce the definitions of consistency, stability and convergence for an SGLM. It will be shown that in SGLMs, stability and consistency together are equivalent to convergence. Also, by introducing a subclass of SGLMs, we construct methods of this subclass up to the maximal order which possess Runge-Kutta stability property and A-stability for implicit ones. [ABSTRACT FROM AUTHOR]

The construction of symmetric and symplectic exponentially-fitted Runge-Kutta methods for the numerical integration of Hamiltonian systems with oscillatory solutions deserves a lot of interest. In previous papers fourth-order and sixth-order symplectic exponentially-fitted integrators of Gauss type, either with fixed or variable nodes, have been derived. In this paper new such integrators of eighth-order are studied and constructed by making use of the six-step procedure of Ixaru and Vanden Berghe (). Numerical experiments for some oscillatory problems are presented. [ABSTRACT FROM AUTHOR]

This paper aims at the application of an optimized two-step hybrid block method for solving boundary value problems with different types of boundary conditions. The proposed approach produces simultaneously approximations at all the grid points after solving an algebraic system of equations. The final approximate solution is obtained through a homotopy-type strategy which is used in order to get starting values for Newton's method. The convergence analysis shows that the proposed method has at least fifth order of convergence. Some numerical experiments such as Bratu's problem, singularly perturbed, and nonlinear system of BVPs are presented to illustrate the better performance of the proposed approach in comparison with other methods available in the recent literature. [ABSTRACT FROM AUTHOR]

For Hamiltonian systems with non-canonical structure matrices, a new family of fourth-order energy-preserving integrators is presented. The integrators take a form of a combination of Runge–Kutta methods and continuous-stage Runge–Kutta methods and feature a set of free parameters that offer greater flexibility and efficiency. Specifically, we demonstrate that by carefully choosing these free parameters, a simplified Newton iteration applied to the integrators of order four can be parallelizable. This results in faster and more efficient integrators compared with existing fourth-order energy-preserving integrators. [ABSTRACT FROM AUTHOR]

Sixty years ago, Butcher (Butcher Math. Soc. 3, 185–201 1963) characterized a natural tabulation of the order conditions for Runge–Kutta methods of order p as an isomorphism from the set of rooted trees having up to p nodes and provided examples of explicit and implicit methods of several orders. Within a few years, Fehlberg (Fehlberg 1968) derived pairs of explicit methods of successive orders that could be implemented efficiently by using the difference of each pair of estimates to control the local error. Unfortunately, Fehlberg's pairs were deficient for quadrature problems. Subsequently, this author (Verner SIAM J. Numer. Anal. 15, 772–790 1978, Verner SIAM J. Numer. Anal. 16, 857–875 1979) derived parametric families of explicit Runge–Kutta pairs of increasing orders 6 to 9 that avoided this problem altogether. These, and most known explicit methods, have been derived by exploiting certain "simplifying conditions" suggested by Butcher (Butcher Math. Soc. 3, 185–201 1963, Butcher 1987 [pp. 194 f.]) that imposed constraints on subsets of the coefficients and thereby simplified the solution of the order conditions for moderate to high-order methods. "RK-Test-21", a MAPLE program developed recently by Butcher (Butcher 2021), was applied to derive known 13-stage pairs of orders 7 and 8. Unexpectedly, results of this application revealed the existence of some previously unknown pairs—i.e., some with the correct orders that satisfied most, but not all, of the previously known simplifying conditions. This present study develops formulas for directly computing exact coefficients of these new pairs together with others lying within this new parametric family of (13,7-8) pairs. While the best of these new pairs falls short of the best of pairs already known, the properties discovered might be utilized to precisely characterize recently reported higher-order methods found using other approaches by Khashin (Khashin 2009) and Zhang (Zhang 2019) and possibly lead to finding other Runge–Kutta and related yet unknown methods. [ABSTRACT FROM AUTHOR]

This paper proposes and investigates the multidimensional extended Runge-Kutta-Nyström (ERKN) methods for the general second-order oscillatory system y″ + My = f( y, y′) where M is a positive semi-definite matrix containing implicitly the frequencies of the problem. The work forms a natural generalization of our previous work on ERKN methods for the special system y″ + My = f( y) (H. Yang et al. Extended RKN-type methods for numerical integration of perturbed oscillators, Comput. Phys. Comm. 180 (2009) 1777-1794 and X. Wu et al., ERKN integrators for systems of oscillatory second-order differential equations, Comput. Phys. Comm. 181 (2010) 1873-1887). The new ERKN methods, with coefficients depending on the frequency matrix M, incorporate the special structure of the equation brought by the term My into both internal stages and updates. In order to derive the order conditions for the ERKN methods, an extended Nyström tree (EN-tree) theory is established. The results of numerical experiments show that the new ERKN methods are more efficient than the general-purpose RK methods and the adapted RKN methods with the same algebraic order in the literature. [ABSTRACT FROM AUTHOR]

In a recent series of papers, the class of energy-conserving Runge-Kutta methods named Hamiltonian BVMs (HBVMs) has been defined and studied. Such methods have been further generalized for the efficient solution of general conservative problems, thus providing the class of Line Integral Methods (LIMs). In this paper we derive a further extension, which we name Enhanced Line Integral Methods (ELIMs), more tailored for Hamiltonian problems, allowing for the conservation of multiple invariants of the continuous dynamical system. The analysis of the methods is fully carried out and some numerical tests are reported, in order to confirm the theoretical achievements. [ABSTRACT FROM AUTHOR]

It is known that the accuracy in estimating a large number of eigenvalues deteriorates when the standard numerical methods are applied, because of the sharp oscillatory behavior of the corresponding eigenfunctions. One method which has proved to be efficient in treating such problems is the Legendre–Gauss Tau method. In this paper we present an exponentially fitted version of this method and we develop practical formulae to correct the estimated eigenvalues. [ABSTRACT FROM AUTHOR]

In this paper, we study in detail the phase properties and stability of numerical methods for general oscillatory second-order initial value problems whose right-hand side functions depend on both the position y and velocity y . In order to analyze comprehensively the numerical stability of integrators for oscillatory systems, we introduce a novel linear test model y ( t) + ? y( t) + µ y ( t)=0 with µ<2 ?. Based on the new model, further discussions and analysis on the phase properties and stability of numerical methods are presented for general oscillatory problems. We give the new definitions of dispersion and dissipation which can be viewed as an essential extension of the traditional ones based on the linear test model y ( t) + ? y( t)=0. The numerical experiments are carried out, and the numerical results showthatthe analysisofphase properties and stability presentedinthispaper ismoresuitableforthenumericalmethodswhentheyareappliedtothe generaloscillatory second-order initial value problem involving both the position and velocity. [ABSTRACT FROM AUTHOR]

In this paper we consider the family of General Linear Methods (GLMs) for the numerical solution of special second order Ordinary Differential Equations (ODEs) of the type y′′ = f( y( t)), with the aim to provide a unifying approach for the analysis of the properties of consistency, zero-stability and convergence. This class of methods properly includes all the classical methods already considered in the literature (e.g. linear multistep methods, Runge-Kutta-Nyström methods, two-step hybrid methods and two-step Runge-Kutta-Nyström methods) as special cases. We deal with formulation of GLMs and present some general results regarding consistency, zero-stability and convergence. The approach we use is the natural extension of the GLMs theory developed for first order ODEs. [ABSTRACT FROM AUTHOR]

This paper describes a package of computer programs for the unified treatment of initial-value problems for systems of ordinary differential equations. The programs implement a numerical method which is efficient for a general class of differential equations. The user may determine the solutions over finite or infinite intervals. The solutions may have singularities at the end-points of the interval for which the solution is sought. Besides giving the initial values and the analytical expression for the differential equations to be solved the user needs to specify the nature of the singularities and give some other analytical information as described in the paper in order to take advantage of the speed and accuracy of the package described here. [ABSTRACT FROM AUTHOR]

Solving high-order or mixed-order boundary value problems by general purpose software often requires the system to be first converted to a larger equivalent first-order system. The cost of solving such problems is generally O( m 3), where m is the dimension of the equivalent first-order system. In this paper, we show how to reduce this cost by exploiting the special structure the “equivalent” first-order system inherits from the original associated mixed-order system. This technique applies to a broad class of boundary value methods. We illustrate the potential benefits by considering in detail a general purpose Runge–Kutta method and a multiple shooting method. [ABSTRACT FROM AUTHOR]

The paper presents a trigonometrically-fitted implicit third derivative Runge-Kutta-Nystöm method (TTRKNM) whose coefficients depend on the frequency and stepsize for periodic initial value problems. The TTRKNM is a pair of methods which is obtained from its continuous version and applied to produce simultaneous approximations to the solution and its first derivative at each point in the interval of interest. A discussion of the stability property of the method is given. Numerical experiments are performed to demonstrate the accuracy and efficiency of the method. [ABSTRACT FROM AUTHOR]

The purpose of this paper is to analyze the algebraic theory of order for the family of general linear Nyström (GLN) methods introduced in D'Ambrosio et al. (Numer. Algorithm 61(2), 331-349, 2012) with the aim to provide a general framework for the representation and analysis of numerical methods solving initial value problems based on second order ordinary differential equations (ODEs). Our investigation is carried out by suitably extending the theory of B-series for second order ODEs to the case of GLN methods, which leads to a general set of order conditions. This allows to recover the order conditions of numerical methods already known in the literature, but also to assess a general approach to study the order conditions of new methods, simply regarding them as GLN methods: the obtained results are indeed applied to both known and new methods for second order ODEs. [ABSTRACT FROM AUTHOR]

In this paper we define an efficient implementation for the family of low-rank energy-conserving Runge-Kutta methods named Hamiltonian Boundary Value Methods (HBVMs), recently defined in the last years. The proposed implementation relies on the particular structure of the Butcher matrix defining such methods, for which we can derive an efficient splitting procedure. The very same procedure turns out to be automatically suited for the efficient implementation of Gauss-Legendre collocation methods, since these methods are a special instance of HBVMs. The linear convergence analysis of the splitting procedure exhibits excellent properties, which are confirmed by a few numerical tests. [ABSTRACT FROM AUTHOR]

For many systems of differential equations modeling problems in science and engineering, there are natural splittings of the right hand side into two parts, one non-stiff or mildly stiff, and the other one stiff. For such systems implicit-explicit (IMEX) integration combines an explicit scheme for the non-stiff part with an implicit scheme for the stiff part. In a recent series of papers two of the authors (Sandu and Zhang) have developed IMEX GLMs, a family of implicit-explicit schemes based on general linear methods. It has been shown that, due to their high stage order, IMEX GLMs require no additional coupling order conditions, and are not marred by order reduction. This work develops a new extrapolation-based approach to construct practical IMEX GLM pairs of high order. We look for methods with large absolute stability region, assuming that the implicit part of the method is A- or L-stable. We provide examples of IMEX GLMs with optimal stability properties. Their application to a two dimensional test problem confirms the theoretical findings. [ABSTRACT FROM AUTHOR]

The Tau method is a numerical technique that consists in constructing polynomial approximate solutions for ordinary differential equations. This method has two approaches: operational and recursive. The former converts the differential problem to a matrix problem and produces approximations in terms of a prescribed orthogonal polynomials basis. In the recursive approach, we construct approximate solutions in terms of a special set of polynomials { Q( t); k = 0, 1, 2...} called canonical polynomials basis. In some cases, the Q's can be obtained explicitly through a recursive formula. But no analogous formulae are reported in the literature for the general cases. In this paper, utilizing the operational Tau method, we develop an algorithm that allows to generate those canonical polynomials iteratively and explicitly. In addition, we demonstrate the capability of the operational Tau method in treating quadratic optimal control problems governed by ordinary differential equations. [ABSTRACT FROM AUTHOR]