Solving Fractional Dynamical System with Freeplay by Combining Memory-Free Approach and Precise Integration Method

The Yuan-Agrawal (YA) memory-free approach is employed to study fractional dynamical systems with freeplay nonlinearities subjected to a harmonic excitation, by combining it with the precise integration method (PIM). By the YA method, the original equations are transformed into a set of first-order piecewise-linear ordinary differential equations (ODEs). These ODEs are further separated as three linear inhomogeneous subsystems, which are solved by PIM together with a predictor-corrector process. Numerical examples show that the results by the presented method agree well with the solutions obtained by the RungeKutta method and a modified fractional predictor-corrector algorithm. More importantly, the presented method has higher computational efficiency.


Introduction
Recently, fractional derivative (FD) has been widely investigated because such mathematical model has merits over integer derivative in describing complex behaviors of some real systems.It has been confirmed that FD can model the realistic phenomena arising in various disciplines such as physics [1], materials science [2], solid mechanics [3], mechanical vibrations [4], biology [5], economics [6], and control theory [7].
Many numerical solution techniques were developed to obtain analytical, semianalytical, and/or numerical solutions of fractional dynamical systems, for example, the finite difference method [8,9], predictor-corrector approach [10,11], operational matrix method [12,13], variational iteration method [14,15], homotopy perturbation method [16,17], Adomian's decomposition method [18,19], to mention a few.Due to the nonlocal character of the fractional derivative, storing the past responses requires a large amount of computer memory.Moreover, a large amount of computational resources is spent on repeat processing of the convolution that describes FD.Accordingly, it is cumbersome to search for even numerical solutions of fractional dynamical systems in a long time duration.
In order to eliminate the drawback of long memory requirement, Yuan and Agrawal proposed a nonclassical approach [20].In their scheme, the fractional differential equation can be converted into a set of first-order ordinary differential equations which can be solved by the Runge-Kutta (RK) method or the trapezoidal integration rule, and so forth.As there are no convolutions in the converted system, this method is called the Yuan-Agrawal (YA) memory-free approach.Later, Singh and Chatterjee [21] extended the memory-free principle introduced in the YA method, which can improve the computation accuracy.
It is worthy of noting that a troublesome problem in a freeplay model lies in determining switching points.Though the switching points can be approximated as the step length is chosen to be refined enough, the computation is very inefficient because refining the time step will lead to exponentially increasing computation cost.Therefore, it is necessary and worthwhile to propose some more efficient approaches to tackle this problem.The precise integration method (PIM) initiated by Zhong and Williams [22] has been widely applied to various problems modeled by ordinary differential equations such as structural dynamics, optimal control, and flexible multibody dynamics problems [23,24].This method is famous for its high accuracy and computation efficiency.We are motivated by its high precision and efficiency to apply this technique to Bagley-Torvik equation [2] with a piecewise linearity.A predictor-corrector technique will be applied to determine the time for switching from one subsystem to another subsystem.A modified fractional predictor-corrector (MFPC) method [25] will be utilized to validate the presented scheme.Numerical examples show that very accurate numerical results can be provided.Moreover, it is much more efficient than the MFPC or the YA method with RK scheme.

Equations of Motions.
Considering the freeplay nonlinearity of a nonlinear spring, the dynamic equation of a single-degree-of-freedom spring-mass-damping system with fractional derivative is described as where , , and  denote the mass, damping coefficient, and stiffness, respectively,   (), 0 <  < 1, is the derivative of order  of the displacement function (), and  sin() is the externally applied force.Figure 1 shows the sketch of the nonlinear stiffness, (), which is a freeplay nonlinearity usually referred to as a bilinear nonlinearity [26]: where  0 ,   ,   , and  are constants.
Under the transformation of a nondimensional time scale  = , (1) can be transformed into the following equation containing the fractional derivative with respect to nondimensional time : By means of gamma function and Laguerre integral formula, the fractional derivative term in (3) can be eliminated, and it becomes where (, ) and  are introduced variables as Taking derivative of ( 5) with respect to , we can get The integral in (4) can be approximated using the Laguerre integral formula [20] as where   and   ( = 1, . . ., ) are the Laguerre weights and node points, respectively.Therefore, ( 4) and ( 7) can be written as a set of first-order ordinary differential equations as with Define X = [, ], ( 1 , ), . . ., (  , )]  , and ( 9) to (10) can be rewritten as the matrix equations where A  , C  ( = 1, 2, 3), and F are constant matrixes and the superscript denotes derivative with respect to nondimensional time .

Precise Integration Method.
In the precise integration method (PIM), there is a simple yet efficient algorithm to compute the exponential matrix.The small time step, Δ, is further split uniformly as  = Δ/2  with  as large positive number (usually as 20 in real practice).By this means, the exponential matrix can be calculated recursively by The matrix T  is introduced to distinguish the higher-order terms from I. Substitution of ( 17) into ( 16) results in The factorization should be iterated  times, that is, for (iter = 0; iter < , iter + +) T  = 2T  + T  ⋅ T  .After these iterations, the exponential matrix for one time step (Δ) can be finally given as T = I + T  .

A Predictor-Corrector Technique to Determine Switching
Points.Note that ( 13) is a piecewise-linear system consisting of three subsystems, and ( 15)-( 18) are based on the fact that the vibration state remains to be located at the same subsystem.Generally, the state will finally leave one subsystem to another as the displacement  passes through   (  =   or   + ), which indicates a switching point, as shown in Figure 2. A predictor-corrector algorithm is applied to accurately find the time when the vibration state is passing one of the switching points.Denote   =   −   and introduce a ratio If  +1 is exactly at the switching point, we have  +1 = 0 and  = 1.It is predicted that the time needed for   to approach the switching point is about Δ.Then the predicted state is further corrected as Repeating ( 19) and ( 20) until  approaches 1 closely enough with a very small tolerance error such as 10 −12 , the vibration state happening at switching points can be detected accurately and efficiently.

A Single-Degree-of-Freedom Dynamical System.
To validate the feasibility and accuracy of the presented algorithm, the parameters in (1) are chosen to be nondimensional values as  = 0.5,  = 1,  = 1,  = 1,  = 1,  0 = 1,   = −1,   = 0,  = 2, and  = 5.And the initial conditions are chosen as  0 = ] 0 = 0.The results obtained by PIM and RK will be compared with those obtained by MFPC.And the results obtained by MFPC with a very small time step, Δ = 0.001, are regarded as the standard solutions.Figure 3 shows the displacements () obtained by PIM and RK with the time step Δ = 0.1 and 0.01, respectively, when the number of Laguerre points is chosen as  = 30.It can be seen that both the PIM and RK results agree well with the standard solutions.
In order to further check the accuracy, the absolute errors of the displacements obtained by PIM, RK, and MFPC with Δ = 0.1 and 0.01 are shown in Figure 4.It shows that the results obtained by PIM and RK can reach almost the same accuracy.The accuracy of both methods should be about 1% for the amplitude of displacement.Figure 5 presents the absolute errors of results using PIM and RK for Δ = 0.01 when the number of Laguerre points is chosen as  = 15, 20, and 30, respectively.It clearly shows that the PIM and RK results converge as the number of Laguerre node points is increasing.
Figure 6 shows the characteristic of the amplitude of the displacement versus the external force amplitude .Obviously, the amplitude is linearly related to the force amplitude.Furthermore, the amplitude-frequency characteristic of the system is also shown in Figure 7.The amplitude reaches its maximum when the frequency  is close to 1.20 and then decreases as the frequency increases.It indicates a resonance as frequency is varying.By the way, both Figures 6 and 7 verify the accuracy and availability of the PIM results.
The computing times for the PIM, RK, and MFPC to calculate the responses of the system in 100 seconds are shown  in Table 1.The computing times for PIM are much less than that for RK and MFPC.Moreover, the PIM computing times increase very slightly but the RK computing times roughly increase linearly with the increasing number of Laguerre node points.

A Two-Degrees-of-Freedom Dynamical System.
In order to further examine the effectiveness, the presented algorithm is employed to solve a two-degrees-of-freedom dynamical system with fractional derivatives as where M, C, and K represent the mass, damping, and stiffness matrixes, respectively, X = [ 1 ,  2 ]  is the displacement matrix, F is the amplitude matrix of the applied harmonic force, and G(X) means the freeplay nonlinearity.While considering the freeplay nonlinearity of one degree of freedom, such as  1 , G(X) can be expressed as    8.The PIM and RK results coincide very well with the MFPC results, respectively.Figure 9 shows the absolute errors of

)Figure 2 :
Figure 2: Illustration of the predictor-corrector algorithm for detecting the vibration state at a switching point.

Figure 5 :Figure 6 :
Figure 5: Comparison of absolute errors of displacements () obtained by PIM and RK for 15, 20, and 30 Laguerre node points.

Figure 7 :
Figure 7: Comparison of the amplitude-frequency curve versus the external force  obtained by the MFPC and PIM.

Table 1 :
Comparison of the computing times.