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

Authors: Q. X. Liu; X. S. He; J. K. Liu; Y. M. Chen; L. C. Huang
Journal: Mathematical Problems in Engineering
 (2016)

---

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


---

## Body

## 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(1)mD2xt+cDαxt+kMx=fsin⁡ωt,where m, c, and k denote the mass, damping coefficient, and stiffness, respectively, Dαx(t), 0<α<1, is the derivative of order α of the displacement function x(t), and fsin⁡ωt is the externally applied force. Figure 1 shows the sketch of the nonlinear stiffness, M(x), which is a freeplay nonlinearity usually referred to as a bilinear nonlinearity [26]: (2)Mx=M0+x-af,x<afM0+Mfx-af,af≤x≤af+δM0+x-af+δMf-1,x>af+δ,where M0, af, Mf, and δ are constants.Figure 1
Sketch of a freeplay nonlinearity.Under the transformation of a nondimensional time scaleτ=ωt, (1) can be transformed into the following equation containing the fractional derivative with respect to nondimensional time τ: (3)mω2D2xτ+cωαDαxτ+kMx=fsin⁡τ.By means of gamma function and Laguerre integral formula, the fractional derivative term in (3) can be eliminated, and it becomes (4)mω2D2xτ+cωαμ∫0∞φy,τdy+kMx=fsin⁡τ,where φ(y,τ) and μ are introduced variables as(5)φy,τ=y2α-1∫0τe-s-τy2Dxsds,(6)μ=2sin⁡παπ.Taking derivative of (5) with respect to τ, we can get (7)Dφy,τ=y2α-1Dxτ-y2α+1φy,τ.The integral in (4) can be approximated using the Laguerre integral formula [20] as (8)∫0∞φy,τdy≈∑i=1nwieyiφyi,τ,where wi and yi(i=1,…,n) 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 (9)Dxτ=ντ,Dντ=1mω2fsin⁡τ-kmω2Mx-cωαμmω2∑i=1nwieyiφyi,τwith (10)Dφyi,τ=yi2α-1ντ-y2φyi,τ,i=1,2,…,n.Define X=x,ν,φy1,τ,…,φyn,τT, and (9) to (10) can be rewritten as the matrix equations (11)X′=A1X+1mω2Fsin⁡τ+C1,x<afA2X+1mω2Fsin⁡τ+C2,af≤x≤af+δA3X+1mω2Fsin⁡τ+C3,x>af+δ,where Ai,Ci(i=1,2,3), and F 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ξ1τ=cos⁡(τ) and ξ2(τ)=sin⁡(τ) satisfying [27](12)ξ′=-ξ2,ξ2′=ξ1.Rewrite X=x,ν,φy1,τ,…,φyn,τ,ξ1,ξ2T, and (11) becomes (13)X′=BiX+Ci,i=1,2,3,where Bi is the coefficient matrix of subsystem i(i=1,2,3) and Ci is a constant matrix.According to the theories of the ordinary differential equations, the analytical solution of (13) can be given as(14)Xτ=exp⁡τ-τ0·Bi·Xτ0+∫τ0τexp⁡τ-s·Bi·Cids=exp⁡τ-τ0·Bi·Xτ0+Bi-1·exp⁡τ-τ0·Bi-I·Ci.The initial conditions for X(0) are given as x(0)=x0, ν(0)=ν0, φ(yi,0)=0(i=1,…,n), ξ1(0)=1, and ξ2(0)=0. Once a time step (Δτ) is chosen, the solutions for a given time series denoted as [0,Δτ,2Δτ,…,NΔτ] can be generated. For example, at the nth time point, the solution is(15)XnΔτ=exp⁡nΔτ·Bi·Xτ0+Bi-1·exp⁡nΔτ·Bi-I·Ci=exp⁡Δτ·Bin·Xτ0+Bi-1·exp⁡Δτ·Bin-I·Ci.Thus, the key to solve (13) is the computation of the exponential matrix for one time step, denoted as T=exp⁡(Δτ·Bi).


### 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 η=Δτ/2n with n as large positive number (usually as 20 in real practice). By this means, the exponential matrix can be calculated recursively by(16)T=exp⁡Δτ·Bi=exp⁡Bi·Δτ2n2n=exp⁡Bi·η2n.Expanding exp⁡(Bi·η) as a series and retaining several lower-order terms, (17)exp⁡Bi·η≈I+Bi·η+Bi·η22+Bi·η36+Bi·η424=I+Ta.The matrix Ta is introduced to distinguish the higher-order terms from I. Substitution of (17) into (16) results in (18)T=I+Ta2n=I+Ta2n-1·I+Ta2n-1.The factorization should be iterated n times, that is, for (iter=0; iter<n, iter++) Ta=2Ta+Ta·Ta. After these iterations, the exponential matrix for one time step (Δτ) can be finally given as T=I+Ta.


### 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 x passes through xf(xf=af or af+δ), which indicates a switching point, as shown in Figure 2.Figure 2
Illustration of the predictor-corrector algorithm for detecting the vibration state at a switching point.A predictor-corrector algorithm is applied to accurately find the time when the vibration state is passing one of the switching points. Denoteσn=xn-xf and introduce a ratio (19)λ=σnσn-σn+1.If xn+1 is exactly at the switching point, we have σn+1=0 and λ=1. It is predicted that the time needed for xn to approach the switching point is about λΔτ. Then the predicted state is further corrected as(20)Xn+1=exp⁡λΔτ·Bi·Xn+Bi-1·exp⁡λΔτ·Bi-I·Ci=Tλ·Xn+Bi-1·Tλ-I·Ci.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.


## 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(1)mD2xt+cDαxt+kMx=fsin⁡ωt,where m, c, and k denote the mass, damping coefficient, and stiffness, respectively, Dαx(t), 0<α<1, is the derivative of order α of the displacement function x(t), and fsin⁡ωt is the externally applied force. Figure 1 shows the sketch of the nonlinear stiffness, M(x), which is a freeplay nonlinearity usually referred to as a bilinear nonlinearity [26]: (2)Mx=M0+x-af,x<afM0+Mfx-af,af≤x≤af+δM0+x-af+δMf-1,x>af+δ,where M0, af, Mf, and δ are constants.Figure 1
Sketch of a freeplay nonlinearity.Under the transformation of a nondimensional time scaleτ=ωt, (1) can be transformed into the following equation containing the fractional derivative with respect to nondimensional time τ: (3)mω2D2xτ+cωαDαxτ+kMx=fsin⁡τ.By means of gamma function and Laguerre integral formula, the fractional derivative term in (3) can be eliminated, and it becomes (4)mω2D2xτ+cωαμ∫0∞φy,τdy+kMx=fsin⁡τ,where φ(y,τ) and μ are introduced variables as(5)φy,τ=y2α-1∫0τe-s-τy2Dxsds,(6)μ=2sin⁡παπ.Taking derivative of (5) with respect to τ, we can get (7)Dφy,τ=y2α-1Dxτ-y2α+1φy,τ.The integral in (4) can be approximated using the Laguerre integral formula [20] as (8)∫0∞φy,τdy≈∑i=1nwieyiφyi,τ,where wi and yi(i=1,…,n) 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 (9)Dxτ=ντ,Dντ=1mω2fsin⁡τ-kmω2Mx-cωαμmω2∑i=1nwieyiφyi,τwith (10)Dφyi,τ=yi2α-1ντ-y2φyi,τ,i=1,2,…,n.Define X=x,ν,φy1,τ,…,φyn,τT, and (9) to (10) can be rewritten as the matrix equations (11)X′=A1X+1mω2Fsin⁡τ+C1,x<afA2X+1mω2Fsin⁡τ+C2,af≤x≤af+δA3X+1mω2Fsin⁡τ+C3,x>af+δ,where Ai,Ci(i=1,2,3), and F 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ξ1τ=cos⁡(τ) and ξ2(τ)=sin⁡(τ) satisfying [27](12)ξ′=-ξ2,ξ2′=ξ1.Rewrite X=x,ν,φy1,τ,…,φyn,τ,ξ1,ξ2T, and (11) becomes (13)X′=BiX+Ci,i=1,2,3,where Bi is the coefficient matrix of subsystem i(i=1,2,3) and Ci is a constant matrix.According to the theories of the ordinary differential equations, the analytical solution of (13) can be given as(14)Xτ=exp⁡τ-τ0·Bi·Xτ0+∫τ0τexp⁡τ-s·Bi·Cids=exp⁡τ-τ0·Bi·Xτ0+Bi-1·exp⁡τ-τ0·Bi-I·Ci.The initial conditions for X(0) are given as x(0)=x0, ν(0)=ν0, φ(yi,0)=0(i=1,…,n), ξ1(0)=1, and ξ2(0)=0. Once a time step (Δτ) is chosen, the solutions for a given time series denoted as [0,Δτ,2Δτ,…,NΔτ] can be generated. For example, at the nth time point, the solution is(15)XnΔτ=exp⁡nΔτ·Bi·Xτ0+Bi-1·exp⁡nΔτ·Bi-I·Ci=exp⁡Δτ·Bin·Xτ0+Bi-1·exp⁡Δτ·Bin-I·Ci.Thus, the key to solve (13) is the computation of the exponential matrix for one time step, denoted as T=exp⁡(Δτ·Bi).


## 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 η=Δτ/2n with n as large positive number (usually as 20 in real practice). By this means, the exponential matrix can be calculated recursively by(16)T=exp⁡Δτ·Bi=exp⁡Bi·Δτ2n2n=exp⁡Bi·η2n.Expanding exp⁡(Bi·η) as a series and retaining several lower-order terms, (17)exp⁡Bi·η≈I+Bi·η+Bi·η22+Bi·η36+Bi·η424=I+Ta.The matrix Ta is introduced to distinguish the higher-order terms from I. Substitution of (17) into (16) results in (18)T=I+Ta2n=I+Ta2n-1·I+Ta2n-1.The factorization should be iterated n times, that is, for (iter=0; iter<n, iter++) Ta=2Ta+Ta·Ta. After these iterations, the exponential matrix for one time step (Δτ) can be finally given as T=I+Ta.


## 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 x passes through xf(xf=af or af+δ), which indicates a switching point, as shown in Figure 2.Figure 2
Illustration of the predictor-corrector algorithm for detecting the vibration state at a switching point.A predictor-corrector algorithm is applied to accurately find the time when the vibration state is passing one of the switching points. Denoteσn=xn-xf and introduce a ratio (19)λ=σnσn-σn+1.If xn+1 is exactly at the switching point, we have σn+1=0 and λ=1. It is predicted that the time needed for xn to approach the switching point is about λΔτ. Then the predicted state is further corrected as(20)Xn+1=exp⁡λΔτ·Bi·Xn+Bi-1·exp⁡λΔτ·Bi-I·Ci=Tλ·Xn+Bi-1·Tλ-I·Ci.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.


## 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 α=0.5, m=1, c=1, k=1, ω=1, M0=1, af=-1, Mf=0, δ=2, and f=5. And the initial conditions are chosen as x0=ν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 x(τ) obtained by PIM and RK with the time step Δτ=0.1 and 0.01, respectively, when the number of Laguerre points is chosen as n=30. It can be seen that both the PIM and RK results agree well with the standard solutions.Figure 3
Displacementsx(τ) obtained by PIM and RK.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 n=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 4
Absolute errors of the displacementsx(τ) obtained by PIM, RK, and MFPC.Figure 5
Comparison of absolute errors of displacementsx(τ) obtained by PIM and RK for 15, 20, and 30 Laguerre node points.Figure6 shows the characteristic of the amplitude of the displacement versus the external force amplitude f. 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.Figure 6
Comparison of the amplitudes versus the external forcef obtained by the MFPC and PIM.Figure 7
Comparison of the amplitude-frequency curve versus the external forceω obtained by the MFPC and PIM.The computing times for the PIM, RK, and MFPC to calculate the responses of the system in 100 seconds are shown in Table1. 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.Table 1
Comparison of the computing times.
Number of Laguerre node points
Nondimensional time step
PIM (second)
RK (second)
MFPC (second)
n
=
15
Δ
τ
=
0.1
0.1217
13.76
4.3314
Δ
τ
=
0.01
1.1055
24.05
381.05
n
=
20
Δ
τ
=
0.1
0.1509
25.24
4.3314
Δ
τ
=
0.01
1.2516
40.77
381.05
n
=
30
Δ
τ
=
0.1
0.2551
59.85
4.3314
Δ
τ
=
0.01
1.8333
72.60
381.05


### 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 as(21)MD2X+CDαX+KX+GX=Fcos⁡ωt,where M, C, and K represent the mass, damping, and stiffness matrixes, respectively, X=x1,x2T 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 x1, G(X) can be expressed as (22)GX=g1x10,g1x1=M0+x1-af,x1<afM0+Mfx1-af,af≤x1≤af+δM0+x1-af+δMf-1,x1>af+δ.In this case, the numerical values are chosen as M=1001, K=3-2-23, C=0.6000.6, F=100, α=0.5, ω=1, M0=1, af=-1, Mf=0.5, and δ=3. Choosing the number of Laguerre points as n=30, the time history curves of displacements x1 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 x1 obtained by MFPC, PIM, and LK with Δτ=0.1 and 0.01, respectively, compared with the results obtained by MFPC with Δτ=0.001. The attained PIM and RK results almost have the same accuracy.Figure 8
Displacementsx1(τ) for system (21) obtained by PIM and RK, respectively.Figure 9
Absolute errors of the displacementsx1(τ) for system (21) obtained by PIM, RK, and MFPC, respectively.Table2 shows the comparison of the computing times of the MFPC, PIM, and RK upon calculating the responses of system (2) over [0,100]. Obviously, the PIM still possesses substantial advantages over the other two methods for highly dimensional systems.Table 2
Comparison of the computing times for solving system (2) over [0, 100].
Number of Laguerre node points
Nondimensional time step
PIM (second)
RK (second)
MFPC (second)
n
=
10
Δ
τ
=
0.1
0.1783
6.486
8.3587
Δ
τ
=
0.01
1.3452
23.11
773.25
n
=
20
Δ
τ
=
0.1
0.5004
27.90
8.3587
Δ
τ
=
0.01
2.6063
32.62
773.25
n
=
30
Δ
τ
=
0.1
0.9130
68.27
8.3587
Δ
τ
=
0.01
4.3939
72.61
773.25


## 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 α=0.5, m=1, c=1, k=1, ω=1, M0=1, af=-1, Mf=0, δ=2, and f=5. And the initial conditions are chosen as x0=ν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 x(τ) obtained by PIM and RK with the time step Δτ=0.1 and 0.01, respectively, when the number of Laguerre points is chosen as n=30. It can be seen that both the PIM and RK results agree well with the standard solutions.Figure 3
Displacementsx(τ) obtained by PIM and RK.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 n=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 4
Absolute errors of the displacementsx(τ) obtained by PIM, RK, and MFPC.Figure 5
Comparison of absolute errors of displacementsx(τ) obtained by PIM and RK for 15, 20, and 30 Laguerre node points.Figure6 shows the characteristic of the amplitude of the displacement versus the external force amplitude f. 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.Figure 6
Comparison of the amplitudes versus the external forcef obtained by the MFPC and PIM.Figure 7
Comparison of the amplitude-frequency curve versus the external forceω obtained by the MFPC and PIM.The computing times for the PIM, RK, and MFPC to calculate the responses of the system in 100 seconds are shown in Table1. 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.Table 1
Comparison of the computing times.
Number of Laguerre node points
Nondimensional time step
PIM (second)
RK (second)
MFPC (second)
n
=
15
Δ
τ
=
0.1
0.1217
13.76
4.3314
Δ
τ
=
0.01
1.1055
24.05
381.05
n
=
20
Δ
τ
=
0.1
0.1509
25.24
4.3314
Δ
τ
=
0.01
1.2516
40.77
381.05
n
=
30
Δ
τ
=
0.1
0.2551
59.85
4.3314
Δ
τ
=
0.01
1.8333
72.60
381.05


## 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 as(21)MD2X+CDαX+KX+GX=Fcos⁡ωt,where M, C, and K represent the mass, damping, and stiffness matrixes, respectively, X=x1,x2T 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 x1, G(X) can be expressed as (22)GX=g1x10,g1x1=M0+x1-af,x1<afM0+Mfx1-af,af≤x1≤af+δM0+x1-af+δMf-1,x1>af+δ.In this case, the numerical values are chosen as M=1001, K=3-2-23, C=0.6000.6, F=100, α=0.5, ω=1, M0=1, af=-1, Mf=0.5, and δ=3. Choosing the number of Laguerre points as n=30, the time history curves of displacements x1 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 x1 obtained by MFPC, PIM, and LK with Δτ=0.1 and 0.01, respectively, compared with the results obtained by MFPC with Δτ=0.001. The attained PIM and RK results almost have the same accuracy.Figure 8
Displacementsx1(τ) for system (21) obtained by PIM and RK, respectively.Figure 9
Absolute errors of the displacementsx1(τ) for system (21) obtained by PIM, RK, and MFPC, respectively.Table2 shows the comparison of the computing times of the MFPC, PIM, and RK upon calculating the responses of system (2) over [0,100]. Obviously, the PIM still possesses substantial advantages over the other two methods for highly dimensional systems.Table 2
Comparison of the computing times for solving system (2) over [0, 100].
Number of Laguerre node points
Nondimensional time step
PIM (second)
RK (second)
MFPC (second)
n
=
10
Δ
τ
=
0.1
0.1783
6.486
8.3587
Δ
τ
=
0.01
1.3452
23.11
773.25
n
=
20
Δ
τ
=
0.1
0.5004
27.90
8.3587
Δ
τ
=
0.01
2.6063
32.62
773.25
n
=
30
Δ
τ
=
0.1
0.9130
68.27
8.3587
Δ
τ
=
0.01
4.3939
72.61
773.25


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



---

*Source: 9430528-2016-04-11.xml*