Improved Short Memory Principle Method for Solving Fractional Damped Vibration Equations

Abstract

In this paper, an improved short memory principle based on the Grünwald–Letnikov definition is proposed and applied in solving fractional vibration differential equations. The improved idea is to adjust the truncation of memory time in short memory principle (SMP) to the truncation of binomial coefficient terms, and the finite coefficients 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 effectiveness.

1. 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/d{x}^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.

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

(1)${}_{t_0}{}^{GL}D_t^{\alpha }f\left(t\right)=\underset{h\to 0}{\lim }\frac{1}{h^{\alpha }}{\displaystyle \sum }_{j=0}^{\left[\left(t-t_0\right)/h\right]}{\left(-1\right)}^j\left(\begin{array}{c}\alpha \\ j\end{array}\right)f\left(t-jh\right)$

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:

(2)${}_{t_0}{}^{GL}D_t^{\alpha }f\left(t\right)=\frac{1}{h^{\alpha }}{\displaystyle \sum }_{j=0}^{\left[\left(t-t_0\right)/h\right]}{w}_{j}f\left(t-jh\right)$

where $w_j={\left(-1\right)}^j \left(\begin{array}{c}\alpha \\ j\end{array}\right)$ is the coefficient of binomial ${\left(1-z\right)}^{\alpha }$:

(3)$w_j=\frac{{{\displaystyle \prod }}_{i=1}^j\left(n-\alpha -1\right)}{j!}, j=1, 2, \cdots$

To avoid the calculation of the Gamma function, the recursive form is adopted:

(4)$w_j=\left(1-\frac{\alpha +1}{j}\right)w_{j-1}, j=1, 2, \cdots$

The initial term is $w_0=\frac{{\left(-1\right)}^0\Gamma \left(\alpha +1\right)}{\Gamma \left(0+1\right)\Gamma \left(\alpha -0+1\right)}=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 $\alpha$. In other words, within $0<\alpha <1$, the memory of the fractional differential decreases with the increase of the derivative order; on the contrary, when $\alpha$ tends to 0, the memory is enhanced. Table 1 and Table 2 provide some binomial coefficients when $\alpha =0.8$ and $\alpha =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₀ = 1 to −1.00 × 10⁻⁴, and the items whose absolute value decreases by approximately 10 times are w₂₀₀, w₈₀₀, w₂₉₃₀, w₁₀₅₃₀ and w₃₇₈₇₀ and attenuate to −1.00 × 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⁻³, and then the items whose absolute value decreases by approximately 10 times are w₄₅₀, w₃₃₈₀, w₂₃₀₄₀, w₁₅₆₉₇₀, w_1069400, w_7225900 and w_45849200. The attenuation to −1.00 × 10⁻¹⁰ requires 45,849,200 items. More than 38 million items from −1.00 × 10⁻⁹ to −1.00 × 10⁻¹⁰, i.e., the truncated items with only 10⁻¹⁰ impact on the current moment will also cause approximately 10⁻² effects. Therefore, the SMP must be improved.

Given the initial conditions, the fractional damped single-degree-of-freedom vibration equation is as follows:

(5)$\left\{\begin{array}{c}m\ddot{x}+c{D}^{\alpha }x+kx=0\\ x_0=0, {\dot{x}}_0=1\end{array}\right.$

where m, c, k and $\alpha$ 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:

(6)$T\left(t\right)={\displaystyle \sum }_{k=0}^q\frac{x_{t0}^{\left(k\right)}}{k!}\left(t-t_0\right)$

where q is the order of the equation. Obviously, $T\left(t\right)=t$ in this example. The equation is constructed as follows:

(7)$x\left(t\right)=z\left(t\right)+T\left(t\right)$

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:

(8)$\left\{\begin{array}{c}m\ddot{z}+c{D}^{\alpha }z+kz=-c{D}^{\alpha }t-kt\\ z_0=0, {\dot{z}}_0=0\end{array}\right.$

According to the derivative definition, step size h is taken to discretise the equation in time and sort it out as follows:

(9)$z_j=\frac{1}{\frac{m}{h^2}+\frac{c}{h^{\alpha }}+k}\left(-\frac{m}{h^2}\left(-2z_{j-1}+z_{j-2}\right)-\frac{c}{h^{\alpha }}{\displaystyle \sum }_{i=1}^j{w}_i{z}_{j-i}-\frac{\mathsf{\Gamma }\left(2\right)}{\mathsf{\Gamma }\left(2-\alpha \right)}{\left(jh\right)}^{\alpha }-k\left(jh\right)\right)$

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 $\alpha$ 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 Figure 1 and Figure 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 Figure 3 and Figure 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 $\alpha =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.

3. 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$\underset{\mathrm{N}\to +\infty }{\lim }\frac{N-\alpha -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\le 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=n{h}_1$, so that the memorable time length becomes $L_2=N_t{\mathrm{h}}_2=n{L}_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=n{h}_2=n^2{h}_1$, and the memorable time length is $L_3=N_t{\mathrm{h}}_3=n{L}_2=n^2{L}_1$.

The calculation formula should be adjusted accordingly. When the number of steps is $k\le N_t$, the expression is

(10)${}_{t_0}{}^{GL}D_t^{\alpha }f\left(k h_1\right) =\frac{1}{h_1{}^{\alpha }}{\displaystyle \sum }_{i=0}^k{w}_{i}f\left(\left(k-i\right)h_1\right)$

When the number of steps is $N_t<k\le n{N}_t$, the expression is

(11)${}_{t_0}{}^{GL}D_t^{\alpha }f\left(k h_1\right) =\frac{1}{h_1{}^{\alpha }}{\displaystyle \sum }_{i=0}^{N_t}{w}_{i}f\left(\left(k-i\right)h_1\right)+\frac{1}{h_2{}^{\alpha }}{\displaystyle \sum }_{i=a}^b{w}_{i}f\left(\left(k-ni\right)h_1\right)$

where a and b are $N_t/n+1$ and $N_t/n+\left[(k-N_t)/n\right]$, respectively, and $\left[·\right]$denotes rounding down. The expression when the number of steps is $n{N}_t<k\le n^2{N}_t$ is

(12)${}_{t_0}{}^{GL}D_t^{\alpha }f\left(k h_1\right) =\frac{1}{h_1{}^{\alpha }}{\displaystyle \sum }_{i=0}^{N_t}{w}_{i}f\left(\left(k-i\right)h_1\right)+\frac{1}{h_2{}^{\alpha }}{\displaystyle \sum }_{i=a}^{N_t}{w}_{i}f\left(\left(k-ni\right)h_1\right)+\frac{1}{h_3{}^{\alpha }}{\displaystyle \sum }_{i=a}^d{w}_{i}f\left(\left(k-n^{2}i\right)h_1\right)$

where d is $N_t/n+\left[(k-n{N}_t)/n^2\right]$.

4. Examples Analysis

4.1. Free Vibration

The single-degree-of-freedom fractional-order-damped free vibration differential equation is selected taking $\alpha =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.

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.

4.2. 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:

(13)$\left\{\begin{array}{c}m\ddot{x}+c{D}^{\alpha }x+kx=\sin \left(t\right)\\ x_0=0, {\dot{x}}_0=0\end{array}\right.$

m, c, k are also taken as 1, $\alpha =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.

4.3. 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:

(14)$\sigma \left(t\right)=\mu \epsilon \left(t\right)+\eta D^{\alpha }\epsilon \left(t\right)$

The equations of motion for each mass are

(15)$$\begin{gathered} m_1{\ddot{x}}_1=\left(x_0-x_1\right)k_e-F_c \\ {m}_2{\ddot{x}}_2=F_c \end{gathered}$$

The stress and strain of viscoelastic material are

(16)$$\begin{gathered} \sigma =\frac{F_c}{A} \\ \epsilon =\left(x_1-x_2\right)/l \end{gathered}$$

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:

(17)$\left\{\begin{array}{l}{m}_1{\ddot{x}}_1+\frac{A\eta }{l}{D}^{\alpha }{x}_1-\frac{A\eta }{l}{D}^{\alpha }{x}_2+\left(k_e+\frac{A\mu }{l}\right)x_1-\frac{A\mu }{l}{x}_2=k_e{x}_0\hfill \\ m_2{\ddot{x}}_2-\frac{A\eta }{l}{D}^{\alpha }{x}_1+\frac{A\eta }{l}{D}^{\alpha }{x}_2-\frac{A\mu }{l}{x}_1+\frac{A\mu }{l}{x}_2=0\hfill \end{array}\right.$

Rewritten into matrix form:

(18)$M\ddot{X}+C{D}^{\mathsf{\alpha }}X+KX=F$

where $M=\left[\begin{array}{cc}{m}_1& 0\\ 0& m_2\end{array}\right]$, $C=\left[\begin{array}{cc}A\eta /l& -A\eta /l\\ -A\eta /l& A\eta /l\end{array}\right]$, $K=\left[\begin{array}{cc}{k}_e+A\mu /l& -A\mu /l\\ -A\mu /l& A\mu /l\end{array}\right]$, $F=\left[\begin{array}{c}{k}_e{x}_0\\ 0\end{array}\right]$ and $X=\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]$.

This article focuses on introducing improved algorithms, so for the convenience of calculation, all parameters except $\alpha$ and $x_0$ are taken as 1. When $\alpha =0.5$ and $x_0=\sin \left(\omega t\right)$ ($\omega$ is $2\mathsf{\pi }$), the step size h = 0.01 is taken under the original definition method, and the displacements of masses 1 and 2 are shown in Figure 15 and Figure 16. The comparison of the two figures shows that the isolated mass m₂ 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₂ 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.

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

Author Contributions

Conceptualization, R.M.; Data curation, R.M.; Formal analysis, R.M.; Funding acquisition, J.H.; Investigation, R.M.; Methodology, R.M.; Project administration, J.H.; Resources, J.H.; Software, R.M.; Supervision, J.H.; Validation, X.Y.; Writing—original draft, R.M.; Writing—review & editing, X.Y. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the National Natural Science Foundation of China [Grant No. 11472133].

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Numerical solutions of original definition method (No truncation) and short memory principle (SMP) method (Nt = 1000, 2000, 3000, 5000) for α=0.8.

Figure 2. Numerical solutions of original definition method (No truncation) and SMP method (Nt = 1000, 2000, 3000, 5000) for α=0.2.

Figure 3. The errors between the SMP numerical solution (Nt = 1000, 2000, 3000, 5000) and the original definition numerical solution (No truncation) for α=0.8.

Figure 4. The errors between the SMP numerical solution (Nt = 1000, 2000, 3000, 5000) and the original definition numerical solution (No truncation) for α=0.2.

Figure 5. Displacement response for different methods (Laplace transform method, original definition method and SMP method).

Figure 6. Errors between numerical solution (original definition method h = 0.01, 0.001 and SMP method h = 0.001 and Nt = 2000) and exact solution (α=0.5).

Figure 7. The schematic of SMP (L1=Nth1).

Figure 8. Schematic diagram of original definition method (k≤Nt).

Figure 9. Schematic diagram of extracting function value by the first enlargement of step size (Nt<k≤nNt).

Figure 10. Schematic diagram of extracting function value by the first enlargement of step size (nNt<k≤n2Nt and k≥n2Nt).

Figure 11. The errors between the two numerical methods (the original defined and the improved SMP) and the Laplace transform method.

Figure 12. Forced vibration response.

Figure 13. Error of the two SMP methods (SMP and improved SMP).

Figure 14. Two-degree-of-freedom model of suspension structure.

Figure 15. Displacement for m1.

Figure 16. Displacement for m2.

Figure 17. Errors of the two SMP methods (SMP and improved SMP) of m2.

Table 1

Partial binomial coefficients for α = 0.8.

j	wj	j	wj	j	wj	j	wj
0	1	10	−0.00298	150	−0.0000212	700	−1.32 × 10−6
1	−0.8	20	−0.000823	200	−0.0000126	800	−1.04 × 10−6
2	−0.08	30	−0.000392	250	−8.44 × 10−6	900	−8.39 × 10−7
3	−0.032	40	−0.000232	300	−6.07 × 10−6	1000	−6.94 × 10−7
4	−0.0176	50	−0.000155	350	−4.60 × 10−6	1100	−5.85 × 10−7
5	−0.0113	60	−0.000111	400	−3.62 × 10-6	1200	−5.00 × 10−7
6	−0.00788	70	−0.0000840	450	−2.92 × 10−6	2930	−1.00 × 10−7
7	−0.00586	80	−0.0000660	500	−2.42 × 10−6	10,530	−1.00 × 10−8
8	−0.00454	90	−0.0000533	550	−2.04 × 10−6	37,870	−1.00 × 10−9
9	−0.00363	100	−0.0000441	600	−1.74 × 10−6	136,140	−1.00 × 10−10

Table 2

Partial binomial coefficients for α = 0.2

j	wj	j	wj	j	wj	j	wj
0	1	10	−0.0110	150	−0.000421	700	−0.0000662
1	−0.2	20	−0.00475	200	−0.000298	800	−0.0000564
2	−0.08	30	−0.00291	250	−0.000228	900	−0.0000490
3	−0.048	40	−0.00206	300	−0.000183	1000	−0.0000432
4	−0.0336	50	−0.00157	350	−0.000152	3380	−1.00 × 10−5
5	−0.0255	60	−0.00126	400	−0.000130	23,040	−1.00 × 10−6
6	−0.0204	70	−0.00105	450	−0.000113	156,970	−1.00 × 10−7
7	−0.0169	80	−0.000895	500	−0.0000992	1,069,400	−1.00 × 10−8
8	−0.0144	90	−0.000777	550	−0.0000884	7,225,900	−1.00 × 10−9
9	−0.0125	100	−0.000685	600	−0.0000797	45,849,200	−1.00 × 10−10