The numerical solutions of IVP’s (Initial Value Problems) associated with perturbed and damped oscillators are a frequent problem in pure and applied sciences such as physics, engineering, celestial mechanics, structural mechanics and the molecular biology among others.
One of the ways that are normally used to solve numerically these problems is through the Runge-Kutta methods and its adaptations    .
The first multistep codes, fixed-step and non-linear, for solving forced oscillators, were published by Stiefel and Bettis  and Bettis  . Stiefel and Scheifele  defined the so-called G-functions. Based on them, these authors constructed a series method that supposes a refinement based on Taylor series. These methods have the important property of integrating harmonic oscillations as a fixed frequency without truncation error. A new method TFSTS (Trigonometrically-Fitted Scheifele Two-Step), based on the Scheifele methods, which verifies this property has been published in .
In , a multistep numerical was developed, which, when contrasted with other methods which use Scheifele functions, has the advantage of incorporating a simple algebraic procedure for calculating the coefficients using recurrence formulas.
During the same period, methods were designed using φ-function series , instead of G-function series, which also integrate the homogeneous problem without truncation error.
In  a family of analytical functions is introduced known as T-functions which are dependent on three parameters. The solution of forced and damped oscillator is expressed as a series of T-functions. Furthermore, the series method of the T-functions is zero, stable and convergent. The series method is extremely precise, however, it has the disadvantage that it needs to be adapted to each specific problem.
In this paper, the transformation of the T-functions series method in the multistep scheme is explained. The calculation of the coefficients of the multistep scheme is made through a recurrent procedure based on the existing relationship between the divided differences and the elemental and complete symmetrical functions . The recurrent calculation of the coefficients permits the multistep scheme to be considered as a VSVO (Variable Step Variable Order) type scheme.
The multistep numerical algorithm, based on T-function series, permits integration of the non-perturbed and damped problem without truncation error.
In this paper, the predictor and corrector methods are designed for the integration of forced and damped oscillators of type , with and , being a parameter of perturbation, usually small, is a known constant frequency and is the damping constant.
In order to cancel the perturbation and carry out the exact integration of the previous IVP, the differential linear operator, , is defined with being a third frequency. In the case of a more complex perturbation, which the operator cannot cancel, simpler expressions of it would be obtained, which would facilitate numerical integration of the IVP.
The resolution of three popular test problems is explained which illustrates the excellent performance of the multistep method. Firstly, we show the resolution of a stiff problem proposed by Lambert  ; secondly, a Duffing oscillator   is treated; and finally, the application of our method to the orbital calculus of an equatorial satellite orbit when the perturbation comes from J2 (zonal harmonics) in B-F variables (Burdet-Ferrándiz variables)     is considered. Its accuracy is compared with the known LSODE (Livermore Solver for Ordinary Differential Equations), MGEAR (C. W. Gear’s Multistep Method) and GEAR (C.W. Gear’s Extrapolation Method) codes implemented in MAPLE (MAThematic PLEasure program).
This paragraph presents a brief description of the T-functions, as well as the series method for the integration of forced and damped oscillators based on them .
2.1. Description of the T-Functions
Let us consider the following IVP:
which formulates an IVP corresponding to a forced and damped oscillator with frequencies, where is a perturbation parameter, usually small.
The solution of (1), obtained for the initial given conditions, is analytic in the interval and the perturbation function
admits in , an absolutely convergent power series development in the following manner:
By applying the differential operator a (1), in order to cancel perturbation, the following superior order equation is obtained:
Taking into account the initial conditions and is obtained:
where is the usual notation of the gradient vector .
A more compact notation is introduced:
The IVP (3), may be described in the following manner:
With the help of Taylor development of given in (2), the IVP (4), may be formulated by the equations:
Definition 1 The solutions of , , are denoted by .
Proposition 1 , with .
Proposition 2 The functions , verify the following recurrence law:
Definition 2 Let be the solutions of with initial conditions , .
Theorem 1 The general solution of with , , , is .
Theorem 2 The general solution for (5) is .
Using the functions , it is possible to extend the functions , in the following manner.
Definition 3 with .
Thus the solution of (5) may be written as:
2.2. T-Functions Series Method
We should consider the IVP (1) and with being the solution, which we consider being analytical in the interval , that is, .
We assume that the approximation to the solution and has been calculated.
The T-functions series method is zero-stable and convergent .
3. Explicit Multistep Method for Perturbed and Damped Oscillators
To implement the multistep explicit method, we will base on the T-functions series method, substituting the derivatives of perturbed function that appear in (6) by expressions in term of divided differences, with some coefficients , elements of a matrix , of those we do not know a recurrence relation  . Once the matrix is known, we will set up a recurrent calculus, through matrix for the explicit method. The study of symmetric polynomials and its relation with the divided differences, will allow us the computation of the matrix .
To make a variable step explicit multistep method of p-steps, we use up to -th order divided difference of a function g, , in the grid values .
The divided differences of perturbed function, g satisfy the identity  :
Denoting by the following matrix, with order:
and choosing , verifies the identity
using a more compact notation is possible to write as
Making a truncation and solving the matrix , it results:
Designating by , it can write:
Substituting (7) in (6) and truncating, it obtains
With , denoting:
we obtain the following formulas, for a multistep explicit method
To construct the method we obtain the elements of the matrix by a recurrent procedure, based on the elemental, , and complete symmetrical functions, . For each the elementary symmetrical function , is the sum of all products of r different variables , being: ,
, and the complete symmetrical function , is defined as the sum of all monomials of total degree r, in the variables that is to say:
where being ,
Particularly and .
Between the divided differences of , that we will denote by and the complete symmetrical polynomial the next relation holds: .
If consider the point , it defines the complete symmetrical functions:
in values with , and . The square matrices of order k, and , are inverse to each other.
As and we can write with . In the particular case, we will get .
The divided differences of one function g satisfy the property:
If , as have order in H, due to the last result, we can write:
Considering and expressing those equalities in a matricial way, we have
and as in the arguments , we can write
As for , if we consider then:
The recurrent form of the matrix is obtained through:
where is a diagonal matrix, such that , with and .
The expressions (9) and (10) allow us to compute the matrix by recurrence, from matrix.
Substituting (11) in (8), we obtain the multistep explicit method.
We take and as the approximate value of the solution and derivative in , with starting values and , respectively.
The formal expression of explicit multistep method is
3.1. Residue Calculation
In the explicit method, the parameter is a common factor of truncation error in each step.
It is assumed that the value calculated for and in is exact, i.e. , , and . The value , is also exact, i.e. , with .
It is necessary to calculate the value of the residues , with .
In the T-functions series methods:
In the multistep explicit method of T-functions:
Similarly for with .
So if the multistep explicit method of T-functions is exact, i.e. it integrates exactly the non-perturbed problem.
3.2. Consistency and Stability
If it is supposed that the function of perturbation is lipschitzian in the variables u and v, and considering the definition of stability given by Stetter, the method is stable . It suffices to apply the theorems 1 and 2 of Grigorieff . According to the theorem 1 of Grigorieff, it is enough to verify that
,satisfies the Lipschitz condition in the domain of its variables and to verify the stability of the method in the non-perturbed case.
The first part follows that the function f is lipschitzian.
Similarly for , and with
To demonstrate the second part, by the theorem 2 of Grigorieff it is enough to test that the set
This arises naturally from the expressions obtained in the construction of the T-functions, since they are continuous functions in and therefore bounded.
For the above results the method is consistent and stable. The method is convergent applying the result given by Stetter , which indicates that consistency and stability imply convergence.
4. Implicit Method for Perturbed and Damped Oscillators
Similarly for the implicit case, the matrix of that we extract the coefficients we will denote as . The matrix , as described in  , is
Designating by , it can write
substituting (12) in (6) and truncating, it results
With , denoting
we obtain the following formulas, for an implicit multistep method
Once the matrix is known, we will set up a recurrent calculus, through matrix for the implicit method. The matrix , as detailed in , is:
The recurrent form of the matrix is then obtained through:
where is a diagonal matrix, such that , with and .
The expressions (14) and (15) allow us to compute the matrix by recurrence, from matrix.
Substituting (16) in (13), we obtain the implicit multistep method.
We take and as the approximate value of the solution and derivative in , with starting values and , respectively.
The formal expression of implicit multistep method is
Analogously to the multistep explicit method, the multistep implicit method of T-functions is convergent.
Predictor Corrector Method for Perturbed and Damped Oscillators
The predictor-corrector method, with variable step size, of p steps for perturbed oscillators is defined like the one which have as predictor the explicit method and as corrector the implicit method, with the previous definitions.
5. An Application for Highly Oscillatory Stiff Problems
In this section, we present an application of the multistep method based in T-functions to the resolution of stiff and highly oscillatory problems. The good behaviour of the method is presented by comparison with other known codes, implemented in the program packaged NUMERIC DSOLVE MAPLE. This program has been developed in an Intel(R) Xeon(r) processor with 12GB of RAM.
This program is useful because permits easy change of the numbers of digits used in calculations and provides convenient graphic capabilities.
5.1. Problem 1
This stiff problem has been selected in order to demonstrate the accuracy of the multistep predictor-corrector method, when the perturbation function is cancelled.
Let’s consider the following IVP stiff problem:
with initial conditions , and solution independent of :
The eigenvalues of the system are −1 and , which enables its degree of stiffness to be regulated. For the case and proposed in  , the stiff problem is expressed as an oscillator, the IVP is:
with solution and derivative , .
Applying the operator , with , in order to cancel the perturbation function, the following IVP is obtained:
Figure 1 and Figure 2, contrast the graph of the decimal logarithm of absolute value of the relative error of the solution and vs t, calculated using the Series method and Predictor-Corrector method with four T-functions, digits
Figure 1. Problem 1. position with four T-functions (Series Method).
Figure 2. Problem 1. position with four T-functions (Multistep Method).
100 and , with the numerical integration codes MGEAR [msteppart], errorper = Float (1, −13), LSODE, , GEAR, errorper = Float (1, −18).
An application is the analysis of an SDOF (Simple Degree of Freedom) subjected to an earthquake excitation.
The importance of an SDOF analysis lies in that it best shows the interdependence between structure and its properties and the duration of earthquake.
5.2. Problem 2
This problem presents one of the Duffing type equations which model electrical and mechanical problems.
The simplest model of a nonlinear vibratory system is that of the free oscillations of a pendulum, which in the case of small oscillations can be studied using a differential Duffing equation, being the basis for the study of many phenomena of nonlinear dynamics.
Micro oscillators are a different option compared to quartz crystals, as they are easier to miniaturize for integration into an electronic circuit, behaving like a Duffing oscillator with an internal resonance.
It is also worth highlighting the application of this chaotic oscillator in the analysis of harmonic SNR (Signals with low signal to Noise Ratio).
The problem proposed by  , which is a linear oscillator with cubic perturbation function, will be considered:
that admits a first integral .
It has been selected this problem to test the behaviour of the multistep predictor-corrector method when the differential operator does not cancel the perturbation. For the case and , applying the operator with to (17), the next IVP is obtained
The T-functions series method has the disadvantage that it needs to be adapted to each specific problem. For the integration of the oscillator (17), it is necessary, prior to the design of the computational algorithm, to establish the recurrence:
The multistep numerical method avoid the adaptation to each specific problem, as occurs in the T-functions series method.
The accuracy of T-functions series method and the multistep method are similar. Figure 3 and Figure 4, contrast the graph of the decimal logarithm of absolute value of the absolute error of the vs t, calculated using the Series method and Predictor-Corrector method with ten T-functions, digits 40 and , with the numerical integration codes MGEAR [msteppart], errorper = Float(1, −15), LSODE, , GEAR, errorper = Float(1, −13).
Figure 3. Problem 2. position with ten T-functions (Series Method).
Figure 4. Problem 2. position with ten T-functions (Multistep Method).
5.3. Problem 3
In Astrodynamics, obtaining precise solutions for the long-term movement of a planetary system or a satellite that orbits a planet is of continuous interest.
In this problem we apply the new multistep method to calculate numerically the position of an Earth artifical satellite   .
To make easier the evaluation of the method we have selected the problem of an equatorial satellite orbit when the perturbation comes from zonal harmonics J2 in B-F variables   expressed using four decupled oscillators with unit frequency.
with initial conditions:
where are the three direction cosines, u is the inverse of the satellite radius using as independent variable the “true anomaly”, c is the angular momentum that it can be considered as a constant, is the reduced mass and e is the eccentricity of the orbit.
This problem is the main problem in the case of an Earth satellite orbiting at low or medium height, since the other acting perturbations are smaller.
It is also integrable and the qualitative behavior of the solution is known and a closed solution is available in terms of special functions, which allows the numerical solution obtained to be contrasted with the exact solution of the problem.
We consider orbits with null e and high eccentricity, .
The first three equations can be integrated exactly, because they are not perturbed. The last equation will be integrated with the predictor-corrector multistep method.
The error has been calculated using the next first integral:
Figure 5. Problem 3. Trajectory error of a circular equatorial J2-perturbed satellite orbit, integrate with seventeen T-functions.
Figure 6. Problem 3. Trajectory error of an equatorial J2-perturbed satellite orbit, integrate with seventeen T-functions and .
value of the absolute error of the vs t, calculated using the Predictor-Corrector method with seventeen T-functions, digits 40 and , with the numerical integration codes MGEAR [msteppart], errorper = Float(1, −15), LSODE, , GEAR, errorper = Float(1, −15).
6. Conclusion and Suggestions
It has been developed a multistep algorithm, based on the T-functions series method, which integrates exactly the non homogeneous problem when it is possible to cancel the perturbation function by a differential operator. If the differential operator does not cancel the perturbation function, it could be obtained a simple expression which would make easier the numerical integration of IVP. The coefficients of this algorithm are calculated by algebraic recurrences, obtained by substituting the derivatives of the perturbation function by divided differences, allowing the design of multistep type codes VSVO. The multistep method is convergent. This method may successfully compete with known integrators, implemented in MAPLE, as shown in the examples, where the gain in precision achieved has been demonstrated for different IVP’s.
 Yang, H., Wu, X. and Xia, J. (2009) Extended RKN-Type Methods for Numerical Integration of Perturbed Oscillators. Computer Physics Communications, 180, 1777-1794.
 You, X., Zhao, J., Yang, H., Fang, Y. and Wu, X. (2014) Order Conditions for RKN Methods Solving General Second-Order Oscillatory Systems. Numerical Algorithms, 66, 147-176.
 You, X., Zhang, Y. and Zhao, J. (2011) Trigonometrically-Fitted Scheifele Two-Step Methods for Perturbed Oscillators. Computer Physics Communications, 182, 1481-1490.
 Reyes, J.A., García-Alonso, F., Ferrándiz, J.M. and Vigo-Aguiar, J. (2007) Numeric Multistep Variable Methods for Perturbed Linear System Integration. Applied Mathematics and Computation, 190, 63-79.
 García-Alonso, F., Reyes, J.A., Ferrándiz, J.M. and Vigo-Aguiar, J. (2009) Accurate Numerical Integration of Perturbed Oscillatory Systems in Two Frequencies. ACM Transactions on Mathematical Software TOMS, 36, Article 21.
 García-Alonso, F., Cortés-Molina, M., Villacampa, Y. and Reyes, J.A. (2016) Accurate Integration of Forced and Damped Oscillators. University Politehnica of Bucharest Scientific Bulletin—Series A, 78, 193-204.
 Martín, P. and Farto, J.M. (1999) Improved Numerical Integration of Perturbed Oscillators via Average. Applied Mathematics and Computation, 99, 129-139.
 Vigo-Aguiar, J. and Ramos, H. (2015) On the Choice of the Frequency in Trigonometrically-Fitted Methods for Periodic Problems. Journal of Computational and Applied Mathematics, 277, 94-105.
 Ferrándiz, J.M. (1988) A General Canonical Transformation Increasing the Number of Variables with Applications to the Two-Body Problem. Celestial Mechanics, 41, 343-357.
 Ferrándiz, J.M., Sansaturio, M.E. and Pojman, J.R. (1992) Increased Accuracy of Computations in the Main Satellite Problem through Linearization Methods. Celestial Mechanics, 53, 347-364.
 Martín, P. and Ferrándiz, J.M. (1997) Multistep Numerical Methods Based on Scheifele G-Functions with Application to Satellite Dynamics. SIAM Journal on Numerical Analysis, 34, 359-375.
 Vigo-Aguiar, J. and Ferrándiz, J.M. (1998) Higher-Order Variable-Step Algorithms Adapted to the Accurate Numerical Integration of Perturbed Oscillators. Computer in Physics, 12, 467-470.
 Van de Vyver, H. (2007) An Adapted Explicit Hybrid Method of Numerov-Type for the Numerical Integration of Perturbed Oscillators. Applied Mathematics and Computation, 186, 1385-1394.