## 4 Numerical Simulation
### Complex Dynamic Response
Figure 20 shows the friction force characteristic of the pneumatic cylinder model in which the solid line denotes the identified curve for extending stroke, whilst the reverse stroke, that is the retracting, is exhibited by the dashed line. It is noteworthy to see that the identified results match well with the simulated ones marked by the square and circle points according to the extending and retracting stroke, respectively. The friction force in the extending stroke is larger than the retracting one. Accordingly, the dynamic equation of the system is rewritten as below:  
\[ \begin{array}{rl}M\ddot{u}+C\dot{u}=&2F_{s}-2F_{s/1}sign(\dot{u})\tan\alpha\\ &+2F_{s/2}sign(\dot{u}\dot{u})\frac{\dot{u}}{\sqrt{(R+r)^{2}-\mu^{2}}}-Mg+F_{e} \cos(\omega t)\end{array}  \tag{42}\]  
where \(F_{s/1}\) and \(F_{s/2}\) are calculated by Equation (39), and \(F_{s}\) is determined by Equation (1). First of all, the amplitude-frequency curves obtained by multi-scale method and fourth-order Runge-Kutta algorithm are compared, as shown in Figure 21. Herein, the frequency of the excited force is swept up slowly from 0 to 10 Rad/s (denoted by the dashed line). In contract to slowly reducing the frequency of 10 Rad/s to 0, the amplitude-frequency curve is exhibited by the dot line. It can be seen that although the down jump point has a difference between two methods due to approximate error and exiting sliding friction between the piston and cylinder, two curves are still in good agreement. Next, the complex behavior of the proposed system will be investigated through the numerical integration for Equation (42) using a fourth-order Runge-Kutta algorithm with  
Figure 21: Multi-scale method compared with numerical integration. Figure 20: Steady-state friction force characteristic. various initial conditions, including velocity and position. The bifurcation diagram of Equation (42) for \(\hat{V}_{t}=7.89\), \(\mu=1.834\), \(\omega\) changed from 1 to 10 rad/s and the same other parameters as in Figure 8. The value of \(\hat{u}\) is determined by using Poincare map with the period \(T=2\pi/\omega\). The simulated result is plotted in Figure 22, showing that if the initial condition is zero for both of the position and velocity, the dynamic response of the system (denoted by filled circles) will be period-2 oscillation as the excited frequency is within 2.2-2.6 rad/s. However, out of this range the solution is period-1 periodic. Meanwhile, the response marked by squares is realized by sweeping the parameter \(\omega\) from 1 to 10 rad/s, for which the initial position and velocity are set at zero only for the first parameter \(\omega\) and the final state of the system will be considered as the initial condition for the next value of \(\omega\). In this case, the oscillation transformed from period-1 to period-3 periodic occurred with \(\omega\) around of 4 and 8 rad/s. Hence, it revealed from this figure that depending on the initial condition, the dynamic response of the system can be bifurcated at A (1.9 rad/s), B (2.8 rad/s), and C (4 rad/s), meaning that in the parameter \(\omega<1.9\), only a period-1 solution exists, but \(1.9<\omega<2.8\) two cases for two period-1 or period-2 solutions may appear. However, if the parameter is within B and C, there are two cases for period-1 periodic oscillation, and after C the dynamic behavior of the system can be a period-1 or period-3 solution. Furthermore, in order to obtain the solution of the dynamic response, the family of the initial conditions named the attractor-basin phase portrait affecting on the dynamic response will be detected in the three cases following. In the first study case, the parameter \(\omega\) is taken account into at value of 2.4 rad/s, while other parameters are simultaneously set as in Figure 22. As shown in Figure 22, it can be seen that the period-1 or period-2 solution may exist depending on the initial condition.