Abstract
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 Runge-Kutta method and a modified fractional predictor-corrector algorithm. More importantly, the presented method has higher computational efficiency.
1. 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.
2. Mathematical Model
2.1. 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, , , is the derivative of order of the displacement function , and 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 , , , 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 asTaking 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 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 , and (9) to (10) can be rewritten as the matrix equations where , and are constant matrixes and the superscript denotes derivative with respect to nondimensional time .
For further simplification and convenience, we expand the dimension by adding two auxiliary variables and satisfying [27]Rewrite , and (11) becomes where is the coefficient matrix of subsystem and is a constant matrix.
According to the theories of the ordinary differential equations, the analytical solution of (13) can be given asThe initial conditions for are given as , , , , and . Once a time step (Δτ) is chosen, the solutions for a given time series denoted as can be generated. For example, at the th time point, the solution isThus, the key to solve (13) is the computation of the exponential matrix for one time step, denoted as
2.2. 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 with as large positive number (usually as 20 in real practice). By this means, the exponential matrix can be calculated recursively byExpanding as a series and retaining several lower-order terms, The matrix is introduced to distinguish the higher-order terms from . Substitution of (17) into (16) results in The factorization should be iterated times, that is, for (; , ) . After these iterations, the exponential matrix for one time step can be finally given as .
2.3. 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 is exactly at the switching point, we have and . It is predicted that the time needed for to approach the switching point is about . Then the predicted state is further corrected asRepeating (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.
3. Numerical Example
3.1. 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 , , , , , , , , , and . And the initial conditions are chosen as . 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, , are regarded as the standard solutions. Figure 3 shows the displacements obtained by PIM and RK with the time step and 0.01, respectively, when the number of Laguerre points is chosen as . 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 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 when the number of Laguerre points is chosen as , 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.
3.2. 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 aswhere , , and represent the mass, damping, and stiffness matrixes, respectively, is the displacement matrix, is the amplitude matrix of the applied harmonic force, and means the freeplay nonlinearity. While considering the freeplay nonlinearity of one degree of freedom, such as , can be expressed as In this case, the numerical values are chosen as , , , , , , , , , and . Choosing the number of Laguerre points as , the time history curves of displacements obtained by MFPC, PIM, and LK with different time steps are shown in Figure 8. The PIM and RK results coincide very well with the MFPC results, respectively. Figure 9 shows the absolute errors of displacements obtained by MFPC, PIM, and LK with and 0.01, respectively, compared with the results obtained by MFPC with . The attained PIM and RK results almost have the same accuracy.
Table 2 shows the comparison of the computing times of the MFPC, PIM, and RK upon calculating the responses of system (2) over . Obviously, the PIM still possesses substantial advantages over the other two methods for highly dimensional systems.
4. Conclusions
This paper has presented an effective algorithm for solving the dynamic responses of fractional dynamical system with a freeplay nonlinearity. According to the YA memory-free method, the system is transformed into a set of piecewise first-order ordinary differential equations without fractional derivative terms. And the transformed system is further separated into three subsystems. The switching points between subsystems are determined accurately by a predictor-corrector algorithm. Compared with the Runge-Kutta method and MFPC, the PIM is able to obtain very accurate solutions much more efficiently.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
This work is supported by the National Natural Science Foundation of China (11572356, 111272361, and 51108472), Guangdong Province Natural Science Foundation (1414050000412), Guangdong Provincial Science & Technology Program of China (2015A020217004), and Fundamental Research Funds for the Central Universities (15lgzd01).