Next Article in Journal
Study of the Vibrational Predissociation of the NeBr2 Complex by Computational Simulation Using the Trajectory Surface Hopping Method
Previous Article in Journal
Stability of Solutions for Parametric Inverse Nonlinear Cost Transportation Problem
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Fernando García-Alonso
,
José Antonio Reyes
and
Mónica Cortés-Molina
*
Department of Applied Mathematics, University of Alicante, 03690 Alicante, Spain
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Mathematics 2020, 8(11), 2028; https://doi.org/10.3390/math8112028
Submission received: 2 October 2020 / Revised: 7 November 2020 / Accepted: 12 November 2020 / Published: 14 November 2020
(This article belongs to the Section Engineering Mathematics)

Abstract

:
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: x ( t ) + A x ( t ) + C x ( t ) = ε F ( x ( t ) , t ) , x ( 0 ) = x 0 , x ( 0 ) = x 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 coefficients 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 first integral. In the first 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.

1. Introduction

Many problems of physics and engineering are modeled by perturbed and damped differential equations systems of the second order: x ( t ) + A x ( t ) + C x ( t ) = ε F ( x ( t ) , t ) , x ( 0 ) = x 0 , x ( 0 ) = x 0 , t [ a , b ] = I .
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.
These 2DOF systems also appear when attempting to reduce the vibration amplitude of a single degree of freedom (SDOF) system, subjected to harmonic excitation when it is connected to a tuned vibration absorber [1,2,3,4,5,6,7].
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].
Different ideas and motivations led several authors to develop special methods for the integration of oscillatory problems. Among the algorithms designed specifically for the numerical resolution of second-order systems are Falkner methods, the Adams extrapolation formulas, Störner–Cowell methods [20,21], and Gauss–Jackson methods [22,23]. Scheifele [24], Stiefel and Scheifele [25], Stiefel and Bettis [26], and Bettis [27,28], followed by Martín and Ferrándiz [29,30,31], Martín and Farto [32], You et al. [33], Vigo–Aguiar and Ferrándiz [34], Reyes et al. [35], and García-Alonso et al. [36], among others, developed numerical methods for the integration of systems of first-order perturbed linear differential equations. These methods can precisely integrate the unperturbed problem. M. Ruggieri and M. P. Speciale [37] carried out a study on perturbed systems of partial differential equations (PDEs) for viscoelastic media, and in some cases, the authors computed approximate solutions by means of the approximate generator of the first-order approximate group of transformations.
This article presents a series method that retains this property and allows for the precise numerical integration of perturbed and damped linear differential systems of the second order. The new method, under certain conditions, will accurately integrate the perturbed and damped problem. The integration is carried out without any previous transformation of the proposed system.
Moreover, this series method is the basis for the construction of multistep explicit, implicit, and predictor–corrector numerical methods.
To this end, a linear operator D + B will be applied to the system of differential equations, where D is a differential operator defined as D = d d t , and B is an appropriate matrix chosen to be able to annul the perturbation function. A family of real functions, { Ψ j } j N , will be defined with values in M ( m , R ) , with which a series method will be constructed. This series method will integrate exactly, with only the first three Ψ functions, an unperturbed system of the third-order equivalent to the given one.
Finally, three problems are solved herein by the method of the Ψ functions series; we contrast the results obtained with the exact solution of the problem and with its first integral. The first example is the integration of a quasi-periodic orbit [38,39,40,41]. The second one is a structural dynamics problem associated with an earthquake [6,7], and the third one studies an equatorial satellite with perturbation J 2 [10,29,42,43].
The good performance of the series method is evidenced by the comparison of the results obtained against other well-known integrators, such as LSODE, ROSENBROCK, GEAR, and DVERK78.

2. Formulation of the Problem

The IVP to be considered is
x + A x + C x = ε F ( x ( t ) , t ) , x ( 0 ) = x 0 , x ( 0 ) = x 0 , t [ a , b ] = I ,
where x : R R m , A , C M ( m , R ) , ε is a perturbation parameter, usually small, and F : R m × R R m .
The components of the vector perturbation field F ( x ( t ) , t ) are F i ( 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 ( t ) = F ( x ( t ) , t ) is an analytical function in I with respect to t. Designating c n and G n ) ( 0 ) , it is possible to write
G ( t ) = F ( x ( t ) , t ) = n = 0 G n ) ( 0 ) n ! t n = n = 0 t n n ! c n .
The IVP (1) can be written using the linear operator derivation D = d d t :
( D 2 + A D + C ) x = ε G ( t ) , x ( 0 ) = x 0 , x ( 0 ) = x 0 , t [ a , b ] = I .
Let us make B M ( m , R ) ; if D + B is applied in (3), it obtains the IVP:
D 3 x + ( A + B ) D 2 x + ( C + B A ) D x + B C x = ( D + B ) ε G ( t ) ,
where the initial values are
x ( 0 ) = x 0 , x ( 0 ) = x 0 , x ( 0 ) = A x 0 C x 0 + ε G ( 0 ) ,
which has the exact same solution as (1) and (3).
By defining L 3 = D 3 + ( A + B ) D 2 + ( C + B A ) D + B C , if A + B = R , C + B A = S , and B C = T , L 3 can be written as L 3 = D 3 + R D 2 + S D + T .
Henceforth, the IVP (4) is
L 3 ( x ) = ( D + B ) ε G ( t ) , x ( 0 ) = x 0 , x ( 0 ) = x 0 , x ( 0 ) = A x 0 C x 0 + ε G ( 0 ) .
The idea that leads to considering (5) is to eliminate the perturbation with the ( D + B ) operator.
From (2) and (5),
L 3 ( x ) = ε ( D + B ) G ( t ) = ε D G ( t ) + ε B G ( t ) = ε n = 1 t n 1 ( n 1 ) ! c n + ε n = 0 t n n ! B c n = ε n = 0 t n n ! c n + 1 + ε n = 0 t n n ! B c n = ε n = 0 t n n ! c n + 1 + B c n .

3. The Ψ Functions

For the construction of the Ψ j functions, the following IVP is considered:
U j + R U j + S U j + T U j = t j j ! I m   with   j = 0 , 1 , , U j ( 0 ) = 0 ̲ , U j ( 0 ) = 0 ̲ m , U j ( 0 ) = 0 ̲ m ,
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:
Ψ j ( t ) + R Ψ j ( t ) + S Ψ j ( t ) + T Ψ j ( t ) = t j 3 ( j 3 ) ! I m , j = 3 , 4 , , Ψ j ( 0 ) = Ψ j ( 0 ) = Ψ j ( 0 ) = 0 ̲ .
Proposition 1.
The Ψ functions verify the following:
Ψ j ( t ) = Ψ j 1 ( t )   w i t h   j 4 .
Proof. 
( Ψ j ) ( t ) + R ( Ψ j ) ( t ) + S ( Ψ j ) ( t ) + T ( Ψ j ) ( t ) = D ( Ψ j ( t ) + R Ψ j ( t ) + S Ψ j ( t ) + T Ψ j ( t ) ) = D t j 3 ( j 3 ) ! I m = t j 4 ( j 4 ) ! I m   with   j = 4 , 5 ,
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:
Ψ j 3 ( t ) + R Ψ j 2 ( t ) + S Ψ j 1 ( t ) + T Ψ j ( t ) = t j 3 ( j 3 ) ! I m , j 6 .
Proof. 
Ψ j ( t ) = Ψ j 1 ( t ) is replaced with j 4 in (8). □
Definition 2.
Ψ 0 ( t ) , Ψ 1 ( t ) , and Ψ 2 ( t ) are the solutions of the following unperturbed system:
U + R U + S U + T U = 0 ̲ ,
with the initial conditions:
U ( 0 ) = I m , U ( 0 ) = 0 ̲ , U ( 0 ) = 0 ̲ ,
U ( 0 ) = 0 ̲ , U ( 0 ) = I m , U ( 0 ) = 0 ̲ ,
U ( 0 ) = 0 ̲ , U ( 0 ) = 0 ̲ , U ( 0 ) = I m ,
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:
V ( t ) = U ( t ) , W ( t ) = U ( t ) = V ( t ) , W ( t ) + R W ( t ) + S V ( t ) + T U ( t ) = 0 ̲ ,
with the following initial conditions:
U ( 0 ) = U 0 = I m , V ( 0 ) = U ( 0 ) = U 0 = 0 ̲ m , W ( 0 ) = U ( 0 ) = V ( 0 ) = 0 ̲ m ,
where U ( t ) , V ( t ) , a n d W ( t ) M ( m , R ) .
(16) can be written as
U ( t ) V ( t ) W ( t ) = 0 ̲ m I m 0 ̲ m 0 ̲ m 0 ̲ m I m T S R U ( t ) V ( t ) W ( t ) , U ( 0 ) V ( 0 ) W ( 0 ) = I m 0 ̲ m 0 ̲ m ,
and the following is noted:
Λ ( t ) = U ( t ) V ( t ) W ( t ) M ( 3 m × m , R ) .
Equation (17) can be expressed as follows:
Λ ( t ) = 0 ̲ m I m 0 ̲ m 0 ̲ m 0 ̲ m I m T S R Λ ( t ) , Λ ( 0 ) = I m 0 ̲ m 0 ̲ m .
By carrying out a restructuring of (19), we can write
U ( t ) 0 ̲ m 0 ̲ m V ( t ) 0 ̲ m 0 ̲ m W ( t ) 0 ̲ m 0 ̲ m = 0 ̲ m I m 0 ̲ m 0 ̲ m 0 ̲ m I m T S R U ( t ) 0 ̲ m 0 ̲ m V ( t ) 0 ̲ m 0 ̲ m W ( t ) 0 ̲ m 0 ̲ m , U ( 0 ) 0 ̲ m 0 ̲ m V ( 0 ) 0 ̲ m 0 ̲ m W ( 0 ) 0 ̲ m 0 ̲ m = I m 0 ̲ m 0 ̲ m 0 ̲ m 0 ̲ m 0 ̲ m 0 ̲ m 0 ̲ m 0 ̲ m ,
and the following is defined:
Ω 0 ( t ) = U ( t ) 0 ̲ m 0 ̲ m V ( t ) 0 ̲ m 0 ̲ m W ( t ) 0 ̲ m 0 ̲ m M ( 3 m × 3 m , R )
and
M = 0 ̲ m I m 0 ̲ m 0 ̲ m 0 ̲ m I m T S R .
The system (20) can be expressed as
Ω 0 ( t ) = M Ω 0 ( t ) , Ω 0 ( 0 ) = I m 0 ̲ m 0 ̲ m 0 ̲ m 0 ̲ m 0 ̲ m 0 ̲ m 0 ̲ m 0 ̲ m ,
whose solution is
Ω 0 ( t ) = e M t I m 0 ̲ m 0 ̲ m 0 ̲ m 0 ̲ m 0 ̲ m 0 ̲ m 0 ̲ m 0 ̲ m .
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:
V ( t ) = U ( t ) , W ( t ) = U ( t ) = V ( t ) , W ( t ) + R W ( t ) + S V ( t ) + T U ( t ) = 0 ̲ ,
with the following initial conditions:
U ( 0 ) = U 0 = 0 ̲ , V ( 0 ) = U ( 0 ) = U 0 = I m , W ( 0 ) = V ( 0 ) = 0 ̲ m ,
where U ( t ) , V ( t ) , a n d W ( t ) M ( m , R ) .
If the matrix
Ω 1 ( t ) = U ( t ) 0 ̲ m 0 ̲ m V ( t ) 0 ̲ m 0 ̲ m W ( t ) 0 ̲ m 0 ̲ m M ( 3 m × 3 m , R )
is built, (23) can be expressed as
Ω 1 ( t ) = M Ω 1 ( t ) , Ω 1 ( 0 ) = 0 ̲ m 0 ̲ m 0 ̲ m I m 0 ̲ m 0 ̲ m 0 ̲ m 0 ̲ m 0 ̲ m ,
whose solution is
Ω 1 ( t ) = e M t 0 ̲ m 0 ̲ m 0 ̲ m I m 0 ̲ m 0 ̲ m 0 ̲ m 0 ̲ m 0 ̲ m .
The submatrix formed by the m + 1 to 2 m rows and the first m columns of Ω 1 ( t ) is denoted by Ψ 1 ( t ) .
Similarly, the calculation of Ψ 2 ( t ) defines
Ω 2 ( t ) = U ( t ) 0 ̲ m 0 ̲ m V ( t ) 0 ̲ m 0 ̲ m W ( t ) 0 ̲ m 0 ̲ m M ( 3 m × 3 m , R )
and considers the following system:
Ω 2 ( t ) = M Ω 2 ( t ) , Ω 2 ( 0 ) = 0 ̲ m 0 ̲ m 0 ̲ m 0 ̲ m 0 ̲ m 0 ̲ m I m 0 ̲ m 0 ̲ m ,
whose solution is
Ω 2 ( t ) = e M t 0 ̲ m 0 ̲ m 0 ̲ m 0 ̲ m 0 ̲ m 0 ̲ m I m 0 ̲ m 0 ̲ m .
The submatrix formed by the 2 m + 1 to 3 m rows and the first m columns of Ω 2 ( t ) is denoted by Ψ 2 ( t ) .
Proposition 3.
Ψ 3 ( t ) = Ψ 2 ( t ) .
Proof. 
( Ψ 3 ) ( t ) + R ( Ψ 3 ) ( t ) + S ( Ψ 3 ) ( t ) + T ( Ψ 3 ) ( t ) = D ( ( Ψ 3 ) ( t ) + R ( Ψ 3 ) ( t ) + S ( Ψ 3 ) ( t ) + T ( Ψ 3 ) ( t ) ) = D ( I m ) = 0 ̲ m .
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):
( Ψ 3 ) ( 0 ) = R ( Ψ 3 ) ( 0 ) S ( Ψ 3 ) ( 0 ) T ( Ψ 3 ) ( 0 ) + I m = 0 ̲ m + 0 ̲ m + 0 ̲ m + I m = I m .
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 ) . □
Proposition 3 extends Proposition 1 to j 3 and Proposition 2 to j 5 .
Theorem 1.
The IVP
L 3 ( x ) = 0 ̲ m , x ( 0 ) = x 0 , x ( 0 ) = x 0 , x ( 0 ) = A x ( 0 ) C x 0 + ε G ( 0 )
has the following solution:
Ψ 0 ( t ) x 0 + Ψ 1 ( t ) x 0 + Ψ 2 ( t ) x 0 .
Proof. 
Taking into account Definition 2 and the linearity of the operator L 3 , the following is obtained:
L 3 ( Ψ 0 ( t ) x 0 + Ψ 1 ( t ) x 0 + Ψ 2 ( t ) x 0 ) = ( D 3 Ψ 0 ( t ) + R D 2 Ψ 0 ( t ) + S D Ψ 0 ( t ) + T Ψ 0 ( t ) ) x 0 + ( D 3 Ψ 1 ( t ) + R D 2 Ψ 1 ( t ) + S D Ψ 1 ( t ) + T Ψ 1 ( t ) ) x 0 + ( D 3 Ψ 2 ( t ) + R D 2 Ψ 2 ( t ) + S D Ψ 2 ( t ) + T Ψ 2 ( t ) ) x 0 = 0 ̲ m .
Definition 2 is enough to demonstrate the following initial conditions:
Ψ 0 ( 0 ) x 0 + Ψ 1 ( 0 ) x 0 + Ψ 2 ( 0 ) x 0 = I m x 0 + 0 ̲ m x 0 + 0 ̲ m x 0 = x 0 , Ψ 0 ( 0 ) x 0 + Ψ 1 ( 0 ) x 0 + Ψ 2 ( 0 ) x 0 = 0 ̲ m x 0 + I m x 0 + 0 ̲ m x 0 = x 0 , Ψ 0 ( 0 ) x 0 + Ψ 1 ( 0 ) x 0 + Ψ 2 ( 0 ) x 0 = 0 ̲ m x 0 + 0 ̲ m x 0 + I m x 0 = x 0 .
 □
Theorem 2.
The solution of the IVP (6) is
x ( t ) = Ψ 0 ( t ) x 0 + Ψ 1 ( t ) x 0 + Ψ 2 ( t ) x 0 + ε n = 3 Ψ n ( t ) c n 2 + B c n 3 .
Proof. 
It is required to prove that x ( t ) verifies the following:
L 3 ( x ( t ) ) = ( D + B ) ε G ( t ) x ( 0 ) = x 0 , x ( 0 ) = x 0 , x ( 0 ) = A x ( 0 ) C x 0 + ε G ( 0 ) = x 0 .
For the operator’s linearity and Theorem 1, it is necessary to have
L 3 ( x ( t ) ) = L 3 Ψ 0 ( t ) x 0 + Ψ 1 ( t ) x 0 + Ψ 2 ( t ) x 0 ) + ε n = 3 Ψ n ( t ) c n 2 + B c n 3 = L 3 Ψ 0 ( t ) x 0 + Ψ 1 ( t ) x 0 + Ψ 2 ( t ) x 0 + ε n = 3 L 3 Ψ n ( t ) c n 2 + B c n 3 = 0 ̲ m + ε n = 3 D 3 Ψ n ( t ) + R D 2 Ψ n ( t ) + S D Ψ n ( t ) + T Ψ n ( t ) c n 2 + B c n 3 = ε n = 3 t n 3 ( n 3 ) ! I m c n 2 + B c n 3 = ε n = 0 t n n ! I m c n + 1 + B c n = ( D + B ) ε G ( t ) ,
and Definitions 1 and 2 are enough to the demonstrate the following initial conditions:
x ( 0 ) = Ψ 0 ( 0 ) x 0 + Ψ 1 ( 0 ) x 0 + Ψ 2 ( 0 ) x 0 + ε n = 3 Ψ n ( 0 ) c n 2 + B c n 3 = I m x 0 + 0 ̲ m x 0 + 0 ̲ m x 0 + 0 ̲ m = I m x 0 = x 0 , x ( 0 ) = Ψ 0 ( 0 ) x 0 + Ψ 1 ( 0 ) x 0 + Ψ 2 ( 0 ) x 0 + ε n = 3 Ψ n ( 0 ) c n 2 + B c n 3 = 0 ̲ m x 0 + I m x 0 + 0 ̲ m x 0 + ε n = 3 0 ̲ m c n 2 + B c n 3 = x 0 , x ( 0 ) = Ψ 0 ( 0 ) x 0 + Ψ 1 ( 0 ) x 0 + Ψ 2 ( 0 ) x 0 + ε n = 3 Ψ n ( 0 ) c n 2 + B c n 3 = 0 ̲ m x 0 + 0 ̲ m x 0 + I m x 0 + ε n = 2 0 ̲ m c n 2 + B c n 3 = x 0 = A x 0 C x 0 + ε G ( 0 ) .
 □

4. 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
x ( t ) = k = 0 t k k ! a k ,
by substituting the McLaurin expansion (34) in (1), we obtain
k = 0 t k k ! a k + 2 + k = 0 t k k ! A a k + 1 + k = 0 t k k ! C a k = ε k = 0 t k k ! c k .
Therefore, recurrences can be achieved between the a k and c k :
a k + 2 + A a k + 1 + C a k = ε c k , k 0 , w i t h a 0 = x 0 , a 1 = x 0 .
By Theorem 2 and (36), the solution x ( t ) can be written as
x ( t ) = Ψ 0 ( t ) a 0 + Ψ 1 ( t ) a 1 + Ψ 2 ( t ) a 2 + ε n = 3 Ψ n ( t ) c n 2 + B c n 3 ,
defining
b 0 = a 0 , b 1 = a 1 , b 2 = a 2 , b k = ε ( c k 2 + B c k 3 ) , with   k 3 ,
and substituting (36) in (38):
b k = ε ( c k 2 + B c k 3 ) = a k + A a k 1 + C a k 2 + B ( a k 1 + A a k 2 + C a k 3 ) = a k + R a k 1 + S a k 2 + T a k 3 , with   k 3 .
The solution can then be expressed as
x ( t ) = Ψ 0 ( t ) b 0 + Ψ 1 ( t ) b 1 + Ψ 2 ( t ) b 2 + n = 3 Ψ n ( t ) b n = n = 0 Ψ n ( t ) b n .
In (40), in addition to a more compact notation than in (37), the coefficients b i of the series depend only on a i .
Once the coefficients 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
x 1 = Ψ 0 ( h ) b 0 + Ψ 1 ( h ) b 1 + Ψ 2 ( h ) b 2 + n = 0 m 3 Ψ n + 3 ( h ) b n + 3 .
Once the approximation of the solution and its derivatives at a point t = n h , which are noted as x n , x n , and x n , are calculated, verifying
L 3 ( x ( t ) ) = ( D + B ) ε G ( t ) , x ( n h ) = x n , x ( n h ) = x n , x ( n h ) = A x n C x n + ε G ( n h ) = x n ,
to obtain an approximation to the solution at a point t = ( n + 1 ) h , a change to the independent variable t = τ + n h is made, so
L 3 ( x ( τ ) ) = ( D + B ) ε G ( τ ) , x ( 0 ) = x n , x ( 0 ) = x n , x ( 0 ) = A x n C x n + ε G ( 0 ) = x n ,
which leads to the initial situation.
Considering
F ( x ( τ ) , τ + n h ) = G ( τ ) = k = 0 τ k k ! c k ,   with   c k = d k G ( 0 ) d τ k = d k G ( n h ) d t k ,
the approximation to the solution at a point t = ( n + 1 ) h is
x n + 1 = Ψ 0 ( h ) b 0 + Ψ 1 ( h ) b 1 + Ψ 2 ( h ) b 2 + n = 0 m 3 Ψ n + 3 ( h ) b n + 3 ,   with   a 0 = x n , a 1 = x n , a k + 2 + A a k + 1 + C a k = c k , k 0 ,   and   b 0 = a 0 , b 1 = a 1 , b 2 = a 2 , b k = a k + R a k 1 + S a k 2 + T 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
x ( t ) = Ψ 0 ( t ) x 0 + Ψ 1 ( t ) x 0 + Ψ 2 ( t ) x 0 + ε n = 3 Ψ n ( t ) ( c n 2 + B c n 3 ) .
Carrying out a truncation of m + 1 Ψ functions, with m 3 ,
x m ( t ) = Ψ 0 ( t ) x 0 + Ψ 1 ( t ) x 0 + Ψ 2 ( t ) x 0 + ε n = 0 m 3 Ψ n + 3 ( t ) ( c n + 1 + B c n ) .
Since
L 3 ( x m ( t ) ) = ε n = 0 m 3 t n n ! c n + 1 + B c n ,
the corresponding residue is given by
R m ( t ) = ( D + B ) ε G ( t ) L 3 ( x m ( t ) ) = ε n = 0 t n n ! c n + 1 + B c n ε n = 0 m 3 t n n ! c n + 1 + B c n = ε n = m 2 t n n ! c n + 1 + B c n .
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.

5. 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.

5.1. 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.
Let us consider the following IVP:
x ( t ) + x ( t ) = 10 3 e i α t , x ( 0 ) = 1 , x ( 0 ) = β i ,   with   | β | 1   and   α , β R .
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:
x 1 ( t ) + x 1 ( t ) = 10 3 cos ( α t ) x 2 ( t ) + x 2 ( t ) = 10 3 sin ( α t ) , x 1 ( 0 ) = 1 , x 1 ( 0 ) = 0 , x 2 ( 0 ) = 0 , x 2 ( 0 ) = β ,
for which the analytical solution is
x 1 ( t ) = 1 ε 1 α 2 cos ( t ) + ε 1 α 2 cos ( α t ) x 2 ( t ) = β ε α 1 α 2 sin ( t ) + ε 1 α 2 sin ( α t ) .
The solution represents the motion on a perturbation of a circular orbit in the complex plane.
We express (51) as
x 1 ( t ) x 2 ( t ) + 1 0 0 1 x 1 ( t ) x 2 ( t ) = ε cos ( α t ) sin ( α t ) , x 1 ( 0 ) = 1 , x 1 ( 0 ) = 0 , x 2 ( 0 ) = 0 , x 2 ( 0 ) = β .
The matrix that annihilates the perturbation function is B = 0 α α 0 . Applying the operator ( D + B ) to the system (53) results in
x 1 ( t ) x 2 ( t ) + 0 α α 0 x 1 ( t ) x 2 ( t ) + x 1 ( t ) x 2 ( t ) + 0 α α 0 x 1 ( t ) x 2 ( t ) = 0 0 , x 1 ( 0 ) x 2 ( 0 ) = 1 0 , x 1 ( 0 ) x 2 ( 0 ) = 0 β , x 1 ( 0 ) x 2 ( 0 ) = ε 1 0 .
The integration with α = 0.1 , β = 0.995 , and ε = 10 3 is achieved with the following algorithm:
a 0 = x 1 ( 0 ) x 2 ( 0 ) = 1 0 , a 1 = x 1 ( 0 ) x 2 ( 0 ) = 0 β , a 2 = A a 1 C a 0 + ε G ( 0 ) = x 1 ( 0 ) x 2 ( 0 ) = ε 1 0
From k = 1 up to n, the following is calculated:
x k = Ψ 0 ( h ) a 0 + Ψ 1 ( h ) a 1 + Ψ 2 ( h ) a 2 , x k = Ψ 0 ( h ) a 0 + Ψ 1 ( h ) a 1 + Ψ 2 ( h ) a 2 , a 0 = x k , a 1 = x k , a 2 = A x k C x k + ε G ( k h ) ,
following k.
Figure 1 shows the graph of the logarithm of the relative error module of the solution x ( t ) , calculated by means of (45) with three Ψ functions, step size 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).

5.2. 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
2 m 0 0 m x 1 x 2 + 3 c c c 2 c x 1 x 2 + 4 k 2 k 2 k 3 k x 1 x 2 = 2 m 0 0 m 1 1 u g ( t ) .
where m is the mass, c is the damping coefficient, and k is the stiffness coefficient.
If 2 m 0 0 m 1 1 u g ( t ) is an harmonic matrix forcing function, i.e.,
2 m 0 0 m 1 1 u g ( t ) = F 0 sin ( ω 0 t ) F 0 sin ( ω 0 t ) ,
where ω 0 is the frequency of the disturbance, then the Equation (56) is
2 m 0 0 m x 1 x 2 + 3 c c c 2 c x 1 x 2 + 4 k 2 k 2 k 3 k x 1 x 2 = F 0 sin ( ω 0 t ) F 0 sin ( ω 0 t ) .
Normalizing the Equation (58), the following is obtained:
x 1 x 2 + 2 m 0 0 m 1 3 c c c 2 c x 1 x 2 + 2 m 0 0 m 1 4 k 2 k 2 k 3 k x 1 x 2 = F 0 sin ( ω 0 t ) 2 m F 0 sin ( ω 0 t ) m .
At the initial moment that the seismic movement occurs, it is logical to assume that the structure is at rest, i.e., x 1 ( 0 ) x 2 ( 0 ) = 0 0 , x 1 ( 0 ) x 2 ( 0 ) = 0 0 , and t [ 0 , T ] .
Following Steffensen’s techniques [44,45], a new variable is added: x 3 = F 0 2 m ω 0 cos ( ω 0 t ) .
This variable allows us to eliminate the perturbation function of the IVP (59), obtaining a new IVP:
x 1 x 2 x 3 + 3 c 2 m c 2 m 0 c m 2 c m 0 0 0 0 x 1 x 2 x 3 + 2 k m k m 0 2 k m 3 k m 0 0 0 0 x 1 x 2 x 3 = F 0 sin ( ω 0 t ) 2 m F 0 sin ( ω 0 t ) m F 0 cos ( ω 0 t ) 2 m , x 1 ( 0 ) x 2 ( 0 ) x 3 ( 0 ) = 0 0 F 0 2 m ω 0 , x 1 ( 0 ) x 2 ( 0 ) x 3 ( 0 ) = 0 0 0 , t [ 0 , T ] .
The matrix B = 0 0 1 0 0 2 ω 0 2 0 0 eliminates the disturbance function. Applying the operator ( D + B ) to the system (60) results in
x 1 x 2 x 3 + 3 c 2 m c 2 m 1 c m 2 c m 2 ω 0 2 0 0 x 1 x 2 x 3 + 2 k m k m 0 2 k m 3 k m 0 3 c ω 0 2 2 m c ω 0 2 m 0 x 1 x 2 x 3 + 0 0 0 0 0 0 2 k ω 0 2 m ω 0 2 k m 0 x 1 x 2 x 3 = 0 0 0 , x 1 ( 0 ) x 2 ( 0 ) x 3 ( 0 ) = 0 0 F 0 2 m ω 0 , x 1 ( 0 ) x 2 ( 0 ) x 3 ( 0 ) = 0 0 0 , x 1 ( 0 ) x 2 ( 0 ) x 3 ( 0 ) = 0 0 F 0 ω 0 2 m , t [ 0 , T ] .
which is solved by the next algorithm, with F 0 = 14 k i p , m = 1.8 k i p s 2 i n , c = 6 π 25 , and k = 16 π 2 5 . The frequency of the disturbance function ω 0 = 4 π 3 r a d s is equal to the natural frequency of vibration of the following structure:
a 0 = x 1 ( 0 ) x 2 ( 0 ) x 3 ( 0 ) = 0 0 F 0 2 m ω 0 , a 1 = x 1 ( 0 ) x 2 ( 0 ) x 3 ( 0 ) = 0 0 0 , a 2 = A a 1 C a 0 + G ( 0 ) = x 1 ( 0 ) x 2 ( 0 ) x 3 ( 0 ) = 0 0 F 0 ω 0 2 m .
From k = 1 up to n, the following is calculated:
x k = Ψ 0 ( h ) a 0 + Ψ 1 ( h ) a 1 + Ψ 2 ( h ) a 2 , x k = Ψ 0 ( h ) a 0 + Ψ 1 ( h ) a 1 + Ψ 2 ( h ) a 2 , a 0 = x k , a 1 = x k , a 2 = A x k C x k + G ( k h ) ,
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 ( 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).

5.3. 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.
x 1 ( t ) + x 1 ( t ) = 0 x 2 ( t ) + x 2 ( t ) = 0 u ( t ) + u ( t ) = μ c 2 + 12 J 2 c 2 u 2 ( t )
with the following initial conditions:
x 1 ( π ) = 1 , x 2 ( π ) = 0 , u ( π ) = μ ( 1 e ) c 2 , x 1 ( π ) = 0 , x 2 ( π ) = 1 , u ( π ) = 0 .
Making the change of variable t π = τ , we can rewrite (63) in the following form:
x 1 ( τ ) x 2 ( τ ) u ( τ ) + x 1 ( τ ) x 2 ( τ ) u ( τ ) = 0 0 μ c 2 + 12 J 2 c 2 u 2 ( τ ) , x 1 ( 0 ) x 2 ( 0 ) u ( 0 ) = 1 0 μ ( 1 e ) c 2 , x 1 ( 0 ) x 2 ( 0 ) u ( 0 ) = 0 1 0 , τ [ 0 , T ] .
Applying the operator ( D + B ) to the system (64), where B is the null matrix results in
x 1 x 2 u + x 1 x 2 u = 0 0 24 J 2 c 2 u u , x 1 ( 0 ) x 2 ( 0 ) u ( 0 ) = 1 0 μ ( 1 e ) c 2 , x 1 ( 0 ) x 2 ( 0 ) u ( 0 ) = 0 1 0 , x 1 ( 0 ) x 2 ( 0 ) u ( 0 ) = 0 0 μ e c 2 + 12 J 2 c 2 μ c 2 ( 1 e ) 2 .
This is integrated by means of the following algorithm:
a 0 = x 1 ( 0 ) x 2 ( 0 ) x 3 ( 0 ) = 1 0 μ ( 1 e ) c 2 , a 1 = x 1 ( 0 ) x 2 ( 0 ) x 3 ( 0 ) = 0 1 0 , a 2 = x 1 ( 0 ) x 2 ( 0 ) x 3 ( 0 ) = A a 1 C a 0 + G ( 0 ) = 0 0 μ e c 2 + 12 J 2 c 2 μ c 2 ( 1 e ) 2
from   j = 1   up   to     n   calculates : from   k = 1   up   to     m   calculates : c k = 0 0 24 J 2 c 2 i = 0 k 1 k 1 i a i , 3 a k i , 3 a k + 2 = c k a k   following   k . x j = Ψ 0 ( h ) a 0 + Ψ 1 ( h ) a 1 + Ψ 2 ( h ) a 2 + l = 3 m c l 2 Ψ l , x j = Ψ 0 ( h ) a 0 + Ψ 1 ( h ) a 1 + Ψ 2 ( h ) a 2 + l = 3 m c l 2 Ψ l 1 , a 0 = x j , a 1 = x j , a 2 = x j + 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:
  • for e = 0 : μ c 2 = 20 21 , J 2 c 2 = 10 21000 ;
  • for e = 0.99 : μ c 2 = 100 20895 , J 2 c 2 = 50 20895000 .
Figure 4 and Figure 5 show the graph of the logarithm of the relative error module of the solution x ( t ) , calculated by means of (45), 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).

6. 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.

Author Contributions

Conceptualization, F.G.-A. and J.A.R.; methodology, F.G.-A. and J.A.R.; software, F.G.-A., J.A.R., and M.C.-M.; validation, F.G.-A., J.A.R., and M.C.-M.; formal analysis, F.G.-A., J.A.R., and M.C.-M.; investigation, F.G.-A., J.A.R., and M.C.-M.; data curation, F.G.-A., and J.A.R.; writing—original draft preparation, F.G.-A., J.A.R., and M.C.-M.; writing—review and editing, F.G.-A., J.A.R., and M.C.-M.; supervision, F.G.-A. and J.A.R.; project administration, F.G.-A.; funding acquisition, J.A.R. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare that there is no conflict of interest.

References

  1. Newmark, N.M. A method of computation for structural dynamics. ASCE J. Eng. Mech. Div. 1959, 85, 67–94. [Google Scholar]
  2. Chung, J.; Hulbert, G.M. A family of single-step Houbolt time integration algorithms for structural dynamics. Comput. Methods Appl. Mech. Eng. 1994, 118, 1–11. [Google Scholar] [CrossRef]
  3. Craig, R.R., Jr.; Kurdila, A.J. Fundamentals of Structural Dynamics; John Wiley and Sons Inc.: Hoboken, NJ, USA, 2006. [Google Scholar]
  4. Chopra, A.K. Dynamics of Structures: Theory and Applications to Earthquake Engineering, 3rd ed.; Prentice-Hall: Upper Saddle River, NJ, USA, 2007. [Google Scholar]
  5. Gholampour, A.; Ghassemieh, M. New implicit method for analysis of problems in nonlinear structural dynamics. Appl. Comput. Mech. 2011, 5, 15–20. [Google Scholar]
  6. Hart, G.C.; Wong, K. Structural Dynamics for Structural Engineers; John Wiley & Sons, Inc.: Hoboken, NJ, USA, 1999. [Google Scholar]
  7. Bangash, M.Y.H. Analyses, Numerical Computations, Codified Methods; Springer: Berlin/Heidelberg, Germany, 2011. [Google Scholar]
  8. Kunstaanheimo, P.; Stiefel, E. Perturbation theory of Kepler motion based on spinor Regularization. J. Reine Angew. Math. 1965, 218, 204–219. [Google Scholar] [CrossRef]
  9. Burdet, C.A. Le mouvement Keplerian et les oscillateurs harmoniques. J. Reine Angew. Math. 1969, 238, 71–78. [Google Scholar] [CrossRef]
  10. Ferrándiz, J.M. A general canonical transformation increasing the number of variables with application to the two-body problem. Celest. Mech. 1988, 41, 343–357. [Google Scholar]
  11. Simos, T.E. A phase-fitted Runge-Kutta-Nyström methods for the numerical solution of initial value problems with oscillating solutions. Comput. Phys. Commun. 2009, 180, 1839–1846. [Google Scholar]
  12. Yang, H.; Wu, X.; You, X.; Fang, Y. Extended RKN-type methods for numerical integration of perturbed oscillators. Comput. Phys. Commun. 2009, 180, 1777–1794. [Google Scholar] [CrossRef]
  13. Yang, H.; Zeng, X.; Wu, X.; Ru, Z. A Simplified Nystrom-Tree Theory for Extended Runge-Kutta-Nystrom Integrators Solving. Comput. Phys. Commun. 2014, 185, 2841–2850. [Google Scholar] [CrossRef]
  14. González, A.B.; Farto, J.M.; López, D.J. Reformulation of the RKGM Methods using Scheifele Expansions. Appl. Math. Lett. 2000, 13, 63–66. [Google Scholar] [CrossRef] [Green Version]
  15. Fang, Y.; Wu, X. A New Pair of Explicit ARKN Methods for the Numerical Integration of General Perturbed Oscillators. Appl. Numer. Math. 2007, 57, 166–175. [Google Scholar] [CrossRef]
  16. You, X.; Zhao, J.; Yang, H.; Fang, Y.; Wu, X. Order Conditions for RKN Methods Solving General Second-Order Oscillatory Systems. Numer. Algorithms 2014, 66, 147–176. [Google Scholar] [CrossRef]
  17. Martín, P.; López, D.J.; García, A. Implementation of Falkner method for problems of the form y” = f(x,y). Appl. Math. Comput. 2000, 109, 183–187. [Google Scholar] [CrossRef]
  18. Franco, J.M. New Methods for Oscillatory Systems Based on ARKN Methods. Appl. Numer. Math. 2006, 56, 1040–1053. [Google Scholar] [CrossRef]
  19. Saravi, M. A procedure for solving some second-order linear ordinary differential equations. Appl. Math. Lett. 2012, 25, 408–411. [Google Scholar] [CrossRef] [Green Version]
  20. Dormand, J.R. Numerical Methods for Differential Equations. A Computational Approach; CRC Press: Boca Ratón, FL, USA, 1996. [Google Scholar]
  21. Henrici, P. Discrete Variable Methods in Ordinary Differential Equations; John Wiley & Sons: New York, NY, USA, 1962. [Google Scholar]
  22. Herrick, S. Astrodynamics. Vol. 2: Orbit Correction, Perturbation Theory Integration; Van Nostrand Reinhold: London, UK, 1972. [Google Scholar]
  23. González, A.B.; Martín, P. A note concerning Gauss-Jackson method. Extr. Math. 1996, 11, 255–260. [Google Scholar]
  24. Scheifele, G. On Numerical Integration of Perturbed Linear Oscillating Systems. ZAMP 1971, 22, 186–210. [Google Scholar] [CrossRef]
  25. Stiefel, E.L.; Scheifele, G. Linear and Regular Celestial Mechanics; Springer: New York, NY, USA, 1971. [Google Scholar]
  26. Stiefel, E.; Bettis, D.G. Stabilization of Cowell’s Method. Numer. Math. 1969, 13, 154–175. [Google Scholar] [CrossRef]
  27. Bettis, D.G. Stabilization of finite difference methods of numerical integration. Celest. Mech. 1970, 2, 282–295. [Google Scholar] [CrossRef]
  28. Bettis, D.G. Numerical integration of products of Fourier and ordinary polinomials. Numer. Math. 1970, 14, 421–434. [Google Scholar] [CrossRef]
  29. Martín, P.; Ferrándiz, J.M. Behaviour of the SMF method for the numerical integration of satellite orbits. Celest. Mech. Dyn. Astron. 1995, 63, 29–40. [Google Scholar] [CrossRef]
  30. Martín, P.; Ferrándiz, J.M. Multistep numerical methods based on Scheifele G-functions with application to satellite dynamics. SIAM J. Numer. Anal. 1997, 34, 359–375. [Google Scholar] [CrossRef]
  31. Martín, P.; Ferrándiz, J.M. Numerical integration of perturbed linear systems. Appl. Math. 1999, 31, 183–189. [Google Scholar] [CrossRef]
  32. Martín, P.; Farto, J.M. Increasing the order of the SMF method for a special type of problem. SIAM J. Numer. Anal. 1998, 35, 773–777. [Google Scholar] [CrossRef]
  33. You, X.; Zhang, Y.; Zhao, J. Trigonometrically-Fitted Scheifele Two-Step. Methods for Perturbed Oscillators. Comput. Phys. Commun. 2011, 182, 1481–1490. [Google Scholar] [CrossRef]
  34. Vigo-Aguiar, J.; Ferrándiz, J.M. Higher-order variable-step algorithms adapted to the accurate numerical integration of perturbed oscillators. Comput. Phys. 1998, 12, 467–470. [Google Scholar] [CrossRef] [Green Version]
  35. Reyes, J.A.; García-Alonso, F.; Ferrándiz, J.M.; Vigo-Aguiar, J. Numeric multistep variable methods for perturbed linear system integration. Appl. Math. Comput. 2007, 190, 63–79. [Google Scholar] [CrossRef]
  36. García-Alonso, F.; Reyes, J.A. A new approach for exact integration of some perturbed stiff linear systems of oscillatory type. Appl. Math. Comput. 2009, 215, 2649–2662. [Google Scholar]
  37. Ruggieri, M.; Speciale, M.P. Approximate Analysis of Nonlinear Dissipative Model. Acta Appl. Math. 2014, 132, 549–559. [Google Scholar] [CrossRef]
  38. Vigo-Aguiar, J.; Ferrándiz, J.M. A general procedure for the adaptation of multistep algorithms to the integration of oscillatory problems. SIAM J. Numer. Anal. 1998, 35, 1684–1708. [Google Scholar] [CrossRef]
  39. Simos, T.E.; Vigo-Aguiar, J. Exponentially fitted sympletic integrator. Phys. Rev. E 2003, 67, 1–7. [Google Scholar] [CrossRef] [PubMed]
  40. Ramos, J.I. Piecewise-linearized methods for initial-value problems with oscillating solutions. Appl. Math. Comput. 2006, 181, 123–146. [Google Scholar] [CrossRef]
  41. Van de Vyver, H. Two-step hybrid methods adapted to the numerical integration of perturbed oscillators. arXiv 2006, arXiv:math/0612637v1. [Google Scholar]
  42. Janin, G. Accurate computation of highly eccentric satellite orbits. Celest. Mech. 1974, 10, 451–456. [Google Scholar] [CrossRef]
  43. Ferrándiz, J.M.; Sansaturio, M.E.; Pojman, J.R. Increased accuracy of computations in the main satellite problem through linearization methods. Celest. Mech. 1992, 53, 347–364. [Google Scholar] [CrossRef]
  44. Steffensen, J.F. On the Differential Equations of Hill in the Theory of the Motion of the Moon. Acta Math. 1955, 93, 169–177. [Google Scholar] [CrossRef]
  45. Steffensen, J.F. On the Differential Equations of Hill in the Theory of the Motion of the Moon (II). Acta Math. 1955, 95, 25–37. [Google Scholar] [CrossRef]
Figure 1. The module’s decimal logarithm of the relative error of the equatorial satellite position.
Figure 1. The module’s decimal logarithm of the relative error of the equatorial satellite position.
Mathematics 08 02028 g001
Figure 2. Two-story frame.
Figure 2. Two-story frame.
Mathematics 08 02028 g002
Figure 3. The module’s decimal logarithm of the relative error of the relative displacement between the mass and the ground.
Figure 3. The module’s decimal logarithm of the relative error of the relative displacement between the mass and the ground.
Mathematics 08 02028 g003
Figure 4. Trajectory error of a circular equatorial J 2 -perturbed satellite orbit, integrated with twenty Ψ functions and e = 0 .
Figure 4. Trajectory error of a circular equatorial J 2 -perturbed satellite orbit, integrated with twenty Ψ functions and e = 0 .
Mathematics 08 02028 g004
Figure 5. Trajectory error of an equatorial J 2 -perturbed satellite orbit, integrated with twenty Ψ functions and e = 0.99 .
Figure 5. Trajectory error of an equatorial J 2 -perturbed satellite orbit, integrated with twenty Ψ functions and e = 0.99 .
Mathematics 08 02028 g005
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

García-Alonso, F.; Reyes, J.A.; Cortés-Molina, M. An Algorithm for the Numerical Integration of Perturbed and Damped Second-Order ODE Systems. Mathematics 2020, 8, 2028. https://doi.org/10.3390/math8112028

AMA Style

García-Alonso F, Reyes JA, Cortés-Molina M. An Algorithm for the Numerical Integration of Perturbed and Damped Second-Order ODE Systems. Mathematics. 2020; 8(11):2028. https://doi.org/10.3390/math8112028

Chicago/Turabian Style

García-Alonso, Fernando, José Antonio Reyes, and Mónica Cortés-Molina. 2020. "An Algorithm for the Numerical Integration of Perturbed and Damped Second-Order ODE Systems" Mathematics 8, no. 11: 2028. https://doi.org/10.3390/math8112028

APA Style

García-Alonso, F., Reyes, J. A., & Cortés-Molina, M. (2020). An Algorithm for the Numerical Integration of Perturbed and Damped Second-Order ODE Systems. Mathematics, 8(11), 2028. https://doi.org/10.3390/math8112028

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop