An Algorithm for the Numerical Integration of Perturbed and Damped Second-Order ODE Systems

: A new method of numerical integration for a perturbed and damped systems of linear second-order differential equations is presented. This new method, under certain conditions, integrates, without truncation error, the IVPs (initial value problems) of the type: xxx (cid:48)(cid:48) ( t ) + A xxx (cid:48) ( t ) + C xxx ( t ) = ε FFF ( xxx ( t ) , t ) , xxx ( 0 ) = xxx 0 , xxx (cid:48) ( 0 ) = xxx (cid:48) 0 , t ∈ [ a , b ] = I , which appear in structural dynamics, astrodynamics, and other fields of physics and engineering. In this article, a succession of real functions is constructed with values in the algebra of m × m matrices. Their properties are studied and we express the solution of the proposed IVP through a serial expansion of the same, whose coefﬁcients are calculated by means of recurrences involving the perturbation function. This expression of the solution is used for the construction of the new numerical method. Three problems are solved by means of the new series method; we contrast the results obtained with the exact solution of the problem and with its ﬁrst integral. In the ﬁrst problem, a quasi-periodic orbit is integrated; in the second, a problem of structural dynamics associated with an earthquake is studied; in the third, an equatorial satellite problem when the perturbation comes from zonal harmonics J 2 is solved. The good behavior of the series method is shown by comparing the results obtained against other integrators.


Introduction
Many problems of physics and engineering are modeled by perturbed and damped differential equations systems of the second order: These systems, known in structural dynamics as two-degree of freedom (2DOF) systems, are used to evaluate the motion of a structure with an attached control system. The dynamic response of structures is interesting for the study of earthquakes.
In celestial mechanics, models such as the classic two-body problem and the satellite problem are reduced to a system of oscillators of this type through the transformations of Kunstaanheimo-Stiefel [8] or Burdet-Ferrándiz [9,10].
The Runge-Kutta type methods [11][12][13][14][15][16][17][18] that usually solve these second-order systems require a large number of evaluations of the perturbation function. This is an inconvenience for its implementation on a computer. Other methods for solving second-order systems can be found in Saravi [19].

Formulation of the Problem
The IVP to be considered is where x x x : R → R m , A, C ∈ M(m, R), ε is a perturbation parameter, usually small, and F F F : R m × R → R m . The components of the vector perturbation field F F F(x x x(t), t) are F i (x x x(t), t) with i = 1, . . . , m, and the field is continuous, with continuous derivatives until a certain order satisfies the conditions for the existence and uniqueness of the solution. This type of system is called a perturbed and damped linear system of the second order. It is considered that G G G(t) = F F F(x x x(t), t) is an analytical function in I with respect to t. Designating c c c n and G G G n) (0), it is possible to write The IVP (1) can be written using the linear operator derivation D = d dt : (D 2 + AD + C) x x x = ε G G G(t), x x x(0) = x x x 0 , x x x (0) = x x x 0 , t ∈ [a, b] = I.
Let us make B ∈ M(m, R); if D + B is applied in (3), it obtains the IVP: where the initial values are which has the exact same solution as (1) and (3).
Henceforth, the IVP (4) is The idea that leads to considering (5) is to eliminate the perturbation with the (D + B) operator. From (2) and (5),

The Ψ Functions
For the construction of the Ψ j functions, the following IVP is considered: where U j ∈ M(m, R), with I m and 0 m being the identity matrix and the zero matrix of this algebra, respectively. Definition 1. The Ψ j functions, which verify Ψ j+3 (t) = U j (t) with j ≥ 0, will be called Ψ functions.
Therefore, the Ψ functions are the Ψ j functions, solutions of the following IVPs: Proposition 1. The Ψ functions verify the following: Of the initial three conditions Ψ j (0) = Ψ j (0) = Ψ j (0) = 0. The one that requires a justification is the third, which is obtained from (8). The other two are evident by the definition.
Since Ψ j (t) and Ψ j−1 (t) verify the same IVP, by the existence and uniqueness theorem of solutions, it can be said that Ψ j (t) = Ψ j−1 (t) with j ≥ 4.

Proposition 2.
The Ψ functions verify the following recurrent relation: Definition 2. Ψ 0 (t), Ψ 1 (t), and Ψ 2 (t) are the solutions of the following unperturbed system: with the initial conditions: respectively.
Ψ i (t), i = 0, 1, 2 can be also calculated as exponential functions of a particular matrix.
In the case Ψ 0 (t), the IVP is considered: with the following initial conditions: where U(t), V(t), and W(t) ∈ M(m, R). (16) can be written as and the following is noted: Equation (17) can be expressed as follows: By carrying out a restructuring of (19), we can write and the following is defined: The system (20) can be expressed as whose solution is We will call Ψ 0 (t) the submatrix formed by the first m + 1 rows and the m columns of Ω(t).
In the case where Ψ 1 (t), the following IVP must be taken into consideration: with the following initial conditions: is built, (23) can be expressed as whose solution is The submatrix formed by the m + 1 to 2m rows and the first m columns of Similarly, the calculation of Ψ 2 (t) defines and considers the following system: whose solution is The submatrix formed by the 2m + 1 to 3m rows and the first m columns of Ω 2 (t) is denoted by Ψ 2 (t).
The initial conditions Ψ 3 (0) and Ψ 3 (0) are 0 m according to the definition of the Ψ functions. The third initial condition is achieved from (8): Since Ψ 3 (t) and Ψ 2 (t) verify the same IVP, by the existence and uniqueness theorem of solutions, we can say that Ψ 3 (t) = Ψ 2 (t).
has the following solution: Proof. Taking into account Definition 2 and the linearity of the operator L 3 , the following is obtained: Definition 2 is enough to demonstrate the following initial conditions: Theorem 2. The solution of the IVP (6) is Proof. It is required to prove that x(t) verifies the following: For the operator's linearity and Theorem 1, it is necessary to have and Definitions 1 and 2 are enough to the demonstrate the following initial conditions:

A New Series Method
In this section, a new numerical series method is shown, based on the Ψ functions, by truncating the solution exposed in Theorem 2.
Through expressing the solution of the IVP (1) in the power series by substituting the McLaurin expansion (34) in (1), we obtain Therefore, recurrences can be achieved between the a a a k and c c c k : a a a k+2 + A a a a k+1 + C a a a k = εc c c k , k ≥ 0, with a a a 0 = x x x 0 , a a a 1 = x x x 0 .
By Theorem 2 and (36), the solution x x x(t) can be written as and substituting (36) in (38): = a a a k + A a a a k−1 + C a a a k−2 + B(a a a k−1 + A a a a k−2 + C a a a k−3 ) = a a a k + R a a a k−1 + S a a a k−2 + T a a a k−3 , with k ≥ 3.
The solution can then be expressed as In (40), in addition to a more compact notation than in (37), the coefficients b b b i of the series depend only on a a a i .
Once the coefficients a a a k from k = 0, . . . , m − 2 are calculated, with step size h, the approximation to the solution of (5) at a point t = h is given by Once the approximation of the solution and its derivatives at a point t = nh, which are noted as x x x n , x x x n , and x x x n , are calculated, verifying to obtain an approximation to the solution at a point t = (n + 1)h, a change to the independent variable t = τ + nh is made, so which leads to the initial situation. Considering the approximation to the solution at a point t = (n + 1)h is a a a k+2 + A a a a k+1 + C a a a k = c c c k , k ≥ 0, and  b b b 0 = a a a 0 , b b b 1 = a a a 1 , b b b 2 = a a a 2 , b b b k = a a a k + R a a a k−1 + S a a a k−2 + T a a a k−3 , k ≥ 3.
The expression (45) provides the Ψ series method that integrates the damped and perturbed system (1).

Residue Calculation
The exact solution of IVP (5) is Carrying out a truncation of m + 1 Ψ functions, with m ≥ 3, Since the corresponding residue is given by Therefore, the perturbation parameter ε is a common factor of the truncation error. If ε = 0, that is, if the perturbation disappears in an arbitrary instant of the independent variable t, the Ψ functions series method integrates the IVP (5) without the truncation error.

Applications for the Ψ Functions Series Method
An implementation of this algorithm is used to approximate the solution of certain IVPs, modeled by second-order systems of differential equations, proposed by various authors.
The results obtained by well-known methods and the Ψ functions series method, versus the analytical solution of the IVP, show the good behavior of the new series method based in Ψ functions.

Problem 1
This problem presents an application of the method to a test problem presented in [25], which can also be found in [26,[38][39][40][41], among others.
The problem is solved as a pair of uncoupled equations. Denoting x(t) = x 1 (t) + x 2 (t) i and by substituting in (50), we get the following second-order system: for which the analytical solution is The solution represents the motion on a perturbation of a circular orbit in the complex plane.
We express (51) as The matrix that annihilates the perturbation function is B = 0 α −α 0 . Applying the operator (D + B) to the system (53) results in The integration with α = 0.1, β = 0.995, and ε = 10 −3 is achieved with the following algorithm: From k = 1 up to n, the following is calculated: following k. Figure 1 shows the graph of the logarithm of the relative error module of the solution x x x(t), calculated by means of (45)

Problem 2
In this problem (2DOF) [6,7], the structure of Figure 2 is solved. This structure is submitted to seismic movement in which only horizontal displacement is considered. x 1 (t) = y 1 (t) − u g (t) and x 2 (t) = y 2 (t) − u g (t) is defined as the relative displacement between the mass and the ground, where y 1 (t) and y 2 (t) are the absolute displacements of the masses, and u g and u g are the absolute ground displacement and the absolute ground velocity, respectively.
The equation of motion for the system is where m is the mass, c is the damping coefficient, and k is the stiffness coefficient.
is an harmonic matrix forcing function, i.e., 2m 0 0 m where ω 0 is the frequency of the disturbance, then the Equation (56) is Normalizing the Equation (58), the following is obtained: At the initial moment that the seismic movement occurs, it is logical to assume that the structure is at rest, i.e., Following Steffensen's techniques [44,45], a new variable is added: x 3 = − F 0 2m ω 0 cos(ω 0 t). This variable allows us to eliminate the perturbation function of the IVP (59), obtaining a new IVP: The matrix B =    0 0 1 0 0 2 −ω 2 0 0 0    eliminates the disturbance function. Applying the operator (D + B) to the system (60) results in which is solved by the next algorithm, with F 0 = 14 kip, m = 1.8 kip s 2 in , c = 6π 25 , and k = 16π 2 5 . The frequency of the disturbance function ω 0 = 4π 3 rad s is equal to the natural frequency of vibration of the following structure: From k = 1 up to n, the following is calculated: following k. Figure 3 shows the graph of the module's decimal logarithm of the relative error of the relative displacement between the mass and the ground x x x(t), calculated by means of (45) with three Ψ functions, a step size of h = 0.1, and 50 digits. This result is compared to the numerical integration algorithms: DEVERK78 with a tolerance of 10 −30 , GEAR with errorper = Float(1,−30), LSODE[backfunc] with a tolerance of 10 −30 , and ROSENBROCK with abserr = Float(1,−30).

Problem 3
In this case, the method is applied to the numerical calculus of the position of an Earth artificial satellite with periodic equatorial orbit that is affected by the perturbation of the zonal harmonic J 2 [30,42]. This problem, expressed by the Burdet-Ferrándiz variables (B-F) [10,43], can be formulated by perturbed and decupled oscillators with unit frequency.
The B-F coordinates are the direction cosines x i and the inverse of the satellite radius using, as an independent variable, the "true anomaly", u. The angular momentum, c, can be considered as a constant, the reduced mass is µ, and the orbit eccentricity is e.
Making the change of variable t − π = τ, we can rewrite (63) in the following form: Applying the operator (D + B) to the system (64), where B is the null matrix results in This is integrated by means of the following algorithm: from j = 1 up to n calculates: from k = 1 up to m calculates: a a a k+2 = c c c k − a a a k following k. x x x j = Ψ 0 (h)a a a 0 + Ψ 1 (h)a a a 1 + Ψ 2 (h)a a a 2 + ∑ m l=3 c c c l−2 Ψ l , x x x j = Ψ 0 (h)a a a 0 + Ψ 1 (h)a a a 1 + Ψ 2 (h)a a a 2 + ∑ m l=3 c c c l−2 Ψ l−1 , a a a 0 = x x x j , a a a 1 = x x x j , a a a 2 = −x x x j + G G G(j h), following j.
A serious drawback occurs in the case of very eccentric orbits since the speed of the satellite in perigee is great, and this circumstance, together with the strong curvature of the trajectory, can cause a loss of precision in integration.
Orbits with zero eccentricity and high eccentricity are considered, respectively:

Conclusions
A family of analytical functions, the Ψ functions, has been defined by studying their properties. Based on these functions and under certain hypotheses, a series method has been constructed that allows, without truncation error, the numerical integration of a wide range of problems modeled by means of systems of second-order differential equations, forced and damped. The Ψ function series method exactly integrates the unperturbed and damped problem when it is possible to cancel the perturbation function by means of a differential operator.
This method has the advantage, compared to other known integrators, of being able to transform a second-order forced and damped problem into an undisturbed problem that integrates in a precise manner, by applying a differential operator.
The integration is carried out without any previous transformation of the proposed system. This series method is the basis for the construction of multistep explicit, implicit, and predictor-corrector numerical methods.
The good behavior of the Ψ functions series method was shown in the following three problems we solved, contrasting the exact solution or the first integral with this method and other known integrators implemented in MAPLE (Deverk78, Gear, Lsode, and Rosenbrock): • The integration of a quasi-periodic orbit-in contrast, the exact solution of the problem is used; • A problem of structural dynamics associated with an earthquake-the analytical solution of the problem was used for comparison; • An equatorial satellite with perturbation J 2 , with high and low eccentricity-the process of contrast was carried out using a first integral.