Improved Short Memory Principle Method for Solving Fractional Damped Vibration Equations

: In this paper, an improved short memory principle based on the Grünwald–Letnikov deﬁnition is proposed and applied in solving fractional vibration di ﬀ erential equations. The improved idea is to adjust the truncation of memory time in short memory principle (SMP) to the truncation of binomial coe ﬃ cient terms, and the ﬁnite coe ﬃ cients are repeatedly applied to the step size gradually enlarged. In this method, a very small initial step size is used to meet the accuracy requirements, and then the step size is gradually enlarged to prolong the memory time and reduce the amount of calculation. Examples of free vibration, forced vibration with a single-degree-of-freedom and a vehicle suspension two-degree-of-freedom vibration reduction model verify the method’s accuracy and e ﬀ ectiveness.


Introduction
The general concept of fractional derivative appears nearly at the same time as integer derivative. At this time, Leibniz expressed the derivative of integer order as d n y/dx n where n as an integer. Then, in 1695, he and L'Hospital discussed the significance of the derivative when n = 1/2 [1]. Fractional calculus did not make good progress in its early development due to the lack of a unified and widely accepted definition. Until the mid-19th century, mathematicians proposed some meaningful definitions of fractional calculus, such as the Grünwald-Letnikov, Riemann-Liouville and Caputo definitions, and then the field of fractional calculus was really established. In recent decades, fractional calculus has been widely used in electrochemical processes [2,3], dielectric polarisation [4], bioengineering [5][6][7], chaos [8], energy supply-demand system [9], viscoelastic mechanics [10,11] and other emerging fields.
Viscoelastic materials are widely used in vibration and noise reduction fields, such as automobiles, airplanes and buildings. For example, the elastomeric lag dampers of helicopters are a single-degree-of-freedom damping mechanism. During operation, viscoelastic materials convert mechanical energy into heat energy to achieve the purpose of vibration reduction. Many experiments have shown that the stress relaxation of many viscoelastic materials is non-exponential and have a historical memory property. The traditional integer-order viscoelastic differential constitutive model cannot accurately describe the mechanical behaviour of viscoelastic materials. However, reasonable results can be obtained by using the fractional constitutive model to describe the stress-strain relationship of viscoelastic materials [12][13][14][15][16][17][18][19]. In the past 20 years, many studies have introduced the fractional differential constitutive model into vibration analysis with viscoelastic damping. On the basis of the fractional Kelvin-Voigt and Maxwell models, the free damped vibration of the oscillator was analysed in reference [20]. [21] investigated the random vibration of a frequency-dependent fractional derivative dynamic system and proposed a standard formula for predicting the random vibration response. The vibration of a rectangular plate with fractional damping and the effect of viscosity on the overall vibration response were studied [22]. The dynamic response of a plane inhomogeneous anisotropic body composed of linear viscoelastic materials was investigated, and a numerical method for solving multi-term fractional differential equations was developed [23]. In reference [24], a variable fractional-order Pasternak-type viscoelastic model is proposed, and the governing equation is established. The results show that the model's fractional order has a great influence on deflection and the bending moment. The response of the Euler-Bernoulli beams with fractional order viscoelasticity under quasi-static and dynamic loads and the random vibration response of beams with fractional damping were studied in [25] and [26]. The nonstationary vibrations of linear viscoelastic single-, two-and multi-degree-of-freedom mechanical systems with a fractional derivative instead of an ordinary derivative were investigated [27]. The application of a fractional-derivative-based isolation model in seismic vibration analysis is discussed in [28]. Nonlinear fractional order systems have also been greatly developed [29][30][31][32][33].
Given the history memory property of the time fractional derivative, a large amount of data is involved in the calculation to obtain more accurate results in the numerical solution process. Reference [12] proposes the short memory principle (SMP) which refers to intercepting the most recent time period and ignoring the time period with less impact and farther away. SMP was used to verify that the fractional derivative of the periodic function is still a periodic function [34]. Through the combination of the SMP and the Grünwald-Letnikov definition, the model reference adaptive control of the fractional order system was studied [35]. According to the classical SMP under the Grünwald-Letnikov definition, Wei et al. proposed and studied several novel SMPs [36]. The traditional SMP may bring large errors, for example, the free vibration cannot return to the equilibrium position. Especially when the fractional order tends to 0, the memory time intercepted by the SMP may be extremely large to ensure calculation accuracy.
This paper proposes an improved SMP based on the Grünwald-Letnikov definition. The improvement aims to adjust the truncation of memory time to the truncation of binomial coefficients. To ensure calculation accuracy, a small step size is initially used. If the number of calculated data points exceeds the number of binomial coefficients, the step size will be enlarged for the limited binomial coefficients to cover a large time area. This improvement aims to obtain more accurate results with as little data as possible. In this paper, a single-degree-of-freedom free vibration example is used to verify the accuracy and reliability of the method.

Memory Effect of Fractional Differential
Unlike the integer order differential, time fractional differential has a historical memory which becomes stronger or weaker with the change of the fractional order. For the fractional differential based on the Grünwald-Letnikov definition, this change is directly reflected in the binomial coefficients. The Grünwald-Letnikov definition is where h is the time step size. If h is sufficiently small, then the limit operation can be ignored and changed into the following form: where w j = (−1) j α j is the coefficient of binomial (1 − z) α : To avoid the calculation of the Gamma function, the recursive form is adopted: The initial term is w 0 = (−1) 0 Γ(α+1) Γ(0+1)Γ(α−0+1) = 1. w j is the weighted coefficient of the function value of the time interval of j steps from the current time, i.e., the binomial coefficient. The coefficient on the left side of the equal sign of the recurrence formula which is less than 1 indicates that this weighting coefficient decreases with the increase of the time interval. The weighting coefficient is also related to the order of derivation. At the same step size, the coefficient decreases with the increase of α. In other words, within 0 < α < 1, the memory of the fractional differential decreases with the increase of the derivative order; on the contrary, when α tends to 0, the memory is enhanced. Tables 1 and 2 provide some binomial coefficients when α = 0.8 and α = 0.2, respectively. The two tables indicate that the absolute value of w j tends to decrease as j increases, and the rate of decrease decelerates. Table 1 shows that it takes only approximately 60 items to decay from w 0 = 1 to −1.00 × 10 −4 , and the items whose absolute value decreases by approximately 10 times are w 200 , w 800 , w 2930 , w 10530 and w 37870 and attenuate to −1.00 × 10 −10 , requiring approximately 136,140 items. Although the speed of attenuation decelerates, the limited items intercepted by the SMP are sufficient to ensure calculation accuracy. In Table 2, 70 items are needed to attenuate from w 0 = 1 to −1.00 × 10 −3 , and then the items whose absolute value decreases by approximately 10 times are w 450 , w 3380 , w 23040 , w 156970 , w 1069400 , w 7225900 and w 45849200 . The attenuation to −1.00 × 10 −10 requires 45,849,200 items. More than 38 million items from −1.00 × 10 −9 to −1.00 × 10 −10 , i.e., the truncated items with only 10 −10 impact on the current moment will also cause approximately 10 −2 effects. Therefore, the SMP must be improved.
Given the initial conditions, the fractional damped single-degree-of-freedom vibration equation is as follows: m ..
where m, c, k and α are the mass, the damping coefficient, the stiffness coefficient and the fractional order, respectively. The Taylor auxiliary function is introduced to change the non-zero initial condition into zero initial condition. The Taylor auxiliary function form is as follows: where q is the order of the equation. Obviously, T(t) = t in this example. The equation is constructed as follows: Then, the initial condition of auxiliary function T(t) is the same as that of x(t), and z(t) is a signal with zero initial conditions. By substituting the construction formula into Equation (5), we obtain the following: m ..
According to the derivative definition, step size h is taken to discretise the equation in time and sort it out as follows: The SMP method is used to change the upper limit of the accumulation item on the right to a fixed value (Nt). By substituting the obtained z(t) into Equation (7), the numerical solution of the original equation is obtained.  To facilitate the calculations, m, c and k are all taken as 1, and α is taken as 0.8 and 0.2. The step size is h = 0.001, and the number of truncation items is Nt = 1000, 2000, 3000, 5000; that is, the memory times of the interception are 1, 2, 3 and 5 s, respectively. The numerical results using the SMP and the original definition are shown in Figures 1 and 2.
Given the relatively small magnitude of error, determining whether to return to the equilibrium position on the response curve is difficult. Therefore, the error between the numerical solution using the SMP and the original definition (i.e., the difference between the numerical solution using the SMP minus the original definition solution) is taken for analysis, as shown in Figures 3 and 4. The figures show the following: (1) as the memory time increases, the calculation error decelerates; (2) none tends towards the equilibrium position; (3) the order of the fractional derivative tends to 0, and the same number of truncated steps will bring greater error. Figure 5 shows the response curves of the numerical and analytical solutions for different cases with α = 0.5. The error results between the numerical and analytical solution (Laplace transform method) are shown in Figure 6. Obviously, (1) if no truncation exists, the response curves will return to the equilibrium position, and the results of small steps will be more accurate, but will increase the amount of calculation; (2) the truncation memory time will reduce the amount of calculation, but will bring greater error and will not tend towards the equilibrium position. Therefore, this paper proposes an improved SMP to reduce the amount of calculation while ensuring the accuracy of calculation as much as possible.

Improved SMP
The calculation formula of the binomial coefficient indicates that when the number of terms is sufficiently large, the change is very small because lim N→+∞ N−α−1 N = 1. Especially when the fractional order approaches 0, a long enough memory time is required to ensure the calculation accuracy. This paper is based on the traditional SMP to improve the calculation accuracy and reduce the amount of calculation. This improvement aims to replace the truncation of time with the truncation of the binomial coefficient in the traditional SMP and gradually increase the step size, so that the binomial coefficient of the finite term can cover a sufficiently long time range. The traditional SMP shown in Figure 7 indicates that the memory time length is L 1 = N t h 1 where N t is the number of binomial coefficient terms intercepted.  Figure 8 shows that when calculating the number of steps (k ≤ N t ), the total memory time is less than or equal to the memory time of the SMP, requiring no step size adjustment. The SMP after the first adjustment and amplification is shown in Figure 9. The step size h 1 for the first enlargement h 2 = nh 1 , so that the memorable time length becomes L 2 = N t h 2 = nL 1 , where n is the step size magnification, and n = 2 in the example figures. Given that L 1 has already been calculated, this part need not be calculated after the step size is enlarged.  Figure 10 shows the entire second step enlargement and a part of the third step enlargement. The second step enlargement uses the same magnification, h 3 = nh 2 = n 2 h 1 , and the memorable time length is L 3 = N t h 3 = nL 2 = n 2 L 1 . The calculation formula should be adjusted accordingly. When the number of steps is k ≤ N t , the expression is When the number of steps is N t < k ≤ nN t , the expression is where a and b are N t /n + 1 and N t /n + [(k − N t )/n], respectively, and [·]denotes rounding down.
The expression when the number of steps is where d is N t /n + (k − nN t )/n 2 .

Free Vibration
The single-degree-of-freedom fractional-order-damped free vibration differential equation is selected taking α = 0.5 conveniently finds the analytical solution. The calculation is carried out in five cases for 20 s: (1) step size h = 0.001, no truncation; (2) step size h = 0.0001, no truncation; (3) step size h = 0.001, the number of coefficient truncation items N t = 1000, step size magnification n = 5; (4) step size h = 0.0001, the number of coefficient truncation terms N t = 2000, the step size magnification n = 15; (5) step size h = 0.0001, the number of coefficient truncation terms N t = 1000, and the step size magnification n = 30. The resulting curve nearly coincides with the exact solution (Laplace transform method) because the error obtained by this method is much smaller than the classical short memory principle. Therefore, the error between the numerical results and the exact solution (i.e., numerical solution minus the exact solution) is compared. The calculation error results are shown in Figure 11. All cases tend towards the equilibrium position, and the two lines with step size of h = 0.001 and h = 0.0001 have no coefficient truncation, indicating that the smaller the step size is, the smaller the error will be. Two error curves have h = 0.001, i.e., the red dashed line without coefficient truncation and the green dotted line with coefficient truncation. Although the curves are nearly coincident, the amount of computation with coefficient truncation is much less than that without coefficient truncation. Although the error of the fourth case is not consistent with the fluctuation position of the error of the second case without coefficient truncation, the error amplitude is roughly the same, also showing the effectiveness of coefficient truncation. In the fifth case, although the initial step size is very small, the effective time of the basic step size is only L 1 = N t h 1 = 0.1, and the step magnification is extremely large, so the error is relatively large. Figure 11. The errors between the two numerical methods (the original defined and the improved SMP) and the Laplace transform method.
When using this method, the number of truncated items and the step size magnification must be coordinated with each other. When the number of truncated items is small, the step size magnification must be reduced, and when the number of truncated items is large, the step size magnification can be appropriately increased. Taking the last numerical result of the second and fourth cases as an example, the comparison of the calculation amounts show that 200,000 items participate in the addition calculation without coefficient truncation, but only approximately 4622 items participate in the addition calculation after the truncation.

Forced Vibration
The improved SMP can make the free vibration return to the equilibrium position. Next, look at the forced vibration. The equation is as follows: m, c, k are also taken as 1, α = 0.5. Equation (13) uses the original fractional definition method, the SMP method (Nt = 3000, 10,000) and the improved SMP method (Nt = 3000, n = 10). The time size is h = 0.0001, and the response curves of 30s are calculated as shown in Figure 12. It can be seen that the numerical solution obtained by the improved SMP method almost coincides with the result of the original definition method, while the result obtained by the SMP method is quite different from them. First of all, the peak value of the curve obtained by the SMP method (Nt = 3000) is smaller than the original definition method; secondly, the SMP method also introduces a visible phase difference, which leads to a larger displacement difference at the same time. The numerical curve obtained by the SMP method (Nt = 10,000) almost completely corrects the phase difference, but its peak value is much higher than the curve obtained by the original definition method. Further, get their error curves (SMP Nt = 10,000 and improved SMP Nt = 3000, n = 10), as shown in Figure 13. The improved SMP adds the time beyond the memory time of the classical SMP by adjusting the step size, and the numerical solution which is very close to the original numerical solution can be obtained. In terms of the amount of calculation, the original definition method has 300,000 items involved in the addition operation in the last step of this example, while the improved SMP has only 8400 items. Although the SMP method (Nt = 10,000) participates in the addition calculation more items than the improved SMP method, it still brings large errors.

Two-Degree-of-Freedom Vibration Reduction Model of Vehicle Suspension
The suspension damping mechanism of the vehicle is simplified to a two-degree-of-freedom model, as shown in Figure 14a. The force analysis is shown in Figure 14b. The constitutive model of viscoelastic material adopts the fractional Kelvin-Voigt model:   .. ..
The stress and strain of viscoelastic material are where A and l are the cross-sectional area and height of viscoelastic material, respectively. Substituting Equations (15) and (16) into (14), the equation of motion of the system can be sorted as follows: ..
Rewritten into matrix form: M ..
where M = m 1 0 0 m 2 , C = Aη/l −Aη/l −Aη/l Aη/l , K = k e + Aµ/l −Aµ/l −Aµ/l Aµ/l , F = k e x 0 0 and This article focuses on introducing improved algorithms, so for the convenience of calculation, all parameters except α and x 0 are taken as 1. When α = 0.5 and x 0 = sin(ωt) (ω is 2π), the step size h = 0.01 is taken under the original definition method, and the displacements of masses 1 and 2 are shown in Figures 15 and 16. The comparison of the two figures shows that the isolated mass m 2 has a better vibration reduction effect. Then use the SMP method (Nt = 1000) and the improved SMP method (Nt = 500, n = 10) to calculate Equation (18). The displacement error of m 2 is shown in Figure 17. Obviously, the improved SMP method effectively reduces the error, and the amount of calculation is also greatly reduced. It can be seen that the improved SMP is also applicable in complex systems.

Conclusions
In this paper, an improved SMP is proposed based on the fractional differential equation of the Grünwald-Letnikov definition, and the single-degree-of-freedom fractional-damped free vibration, forced vibration differential equations and vehicle suspension two-degree-of-freedom vibration reduction model are used as examples to verify its validity and reliability. If the traditional SMP is adopted, than the free vibration will not be attributed to the equilibrium position, and the improved SMP eliminates this error by adjusting the step size and taking points to participate in the calculation in the full time domain. Since free vibration has the characteristics of amplitude attenuation, it can be clearly judged whether to return to the equilibrium position in error analysis. The long-time forced vibration example shows the obvious superiority of the algorithm, which uses less calculation and greatly reduces the error. An example of a two-degree-of-freedom vibration reduction model of a vehicle suspension shows the usability of the improved SMP in the analysis of complex structures. The step size is gradually enlarged, so that the amount of calculation is greatly reduced, and the calculation cost is saved while maintaining the accuracy. To ensure the accuracy and reduce the amount of calculation at the same time, is the purpose of the improved method. It is also the purpose of improving the SMP, using basic small steps to ensure calculation accuracy and gradually enlarging the step size to reduce the amount of calculation. In summary, the improved numerical method based on the SMP is effective and feasible.