# Initial-Sensitive Dynamical Behaviors of a Class of Geometrically Nonlinear Oscillators

Authors: Bo Qin; Huilin Shang; Huimin Jiang
Journal: Shock and Vibration
 (2022)

---

## Abstract

The vibrating system of a class of linkage-slider structure is considered, and its initial-sensitive dynamical behaviors such as safe jump, locking instability, and chaos are studied. First, static bifurcation of the dynamical system is discussed. Then, via analyzing the effect of the external excitation on the periodic solutions under primary resonance, it is found that the change of the excitation frequency may lead to bistability and safe jump. Furthermore, it follows from the investigation of the heteroclinic bifurcation that the increase of the external excitation amplitude may lead to locking instability, chaos, and static locking. The results have some potential values in the design of geometrically nonlinear oscillators.


---

## Body

## 1. Introduction

Dynamics of nonlinear oscillators with quadratic or cubic terms has been investigated extensively in the literature. There are also some oscillators whose nonlinearities are irrational terms due to the changes in their geometric configuration [1, 2]. Such oscillatory systems have received increasing attentions for being commonly used in the design and applications of machinery such as energy harvesters, turbine structures, and industrial robots [3–5]. In most of the cases, it is often quite difficult to get the exact periodic solutions; thus, the pursuit of approximate periodic solutions to these nonlinear oscillators has generated many analytical or semi-analytical methods, for instance, harmonic balance method [6] and parameter-expansion method [7]. Lai and Xiang [8] presented a generalized Senator–Bapat perturbation approach to obtain accurate solutions of an oscillating system with an irrational restoring force. Yildirim et al. [9] applied the Hamiltonian approach to get the analytical approximate solution of the nonlinear oscillators with rational and irrational elastic terms. Li et al. [10] considered a strongly irrational oscillatory system with a viscous damping as well as an external harmonic excitation and employed a four-dimensional average method to obtain the perturbed primary responses. Razzak [11] proposed a new noble modified method to obtain the approximate solution of strongly nonlinear oscillator systems with a rational force and an irrational one. Qin and Shang [12] applied the average method directly to get the approximate periodic solutions semi-analytically in the vibrating system of a linkage slider under an irrational tensile force. Amer et al. [13] used the Poincare method of small parameter to achieve the asymptotic solutions of a three-dimensional system of a gyrostat for the case of irrational frequencies. Kenmogne et al. [14] considered a discontinuous elastic coupling electromechanical system with strong irrational nonlinearity and studied the oscillatory bursting of its electric circuit and the impulse bursting of its mechanical structure via numerical approaches.Apart from the periodic oscillation, complex dynamics of the systems with geometric nonlinearities such as multistability [15, 16], period-n motion [17], and chaos [18] also attracts great attentions during these decades. Among these complicated behaviors, multistability is characterized by a high degree of complexity, as only one of coexisting attractors can guarantee the desired behavior. Cao et al. [19] investigated the global bifurcations and multiple bucklings of a nonlinear oscillator with a pair of strong irrational nonlinear restoring forces theoretically. Santhosh et al. [20] investigated multistability of a harmonically excited SD (smooth and discontinuous) oscillator via numerical integration methods and found that the fractal domains of attraction of the coexisting attractors led to the phenomenon safe jump. Hao and Cao [21] presented a single-degree-of-freedom geometrically nonlinear oscillator with stable quasi-zero stiffness and carried out numerical simulations to detect chaos and coexisting period-3 attractors. Yue et al. [22] studied stochastic bifurcations of the SD oscillator with bounded noises by the generalized cell mapping method and described the bifurcations as a sudden change in shape and size of random attractor with the continuous variation of system parameters. From a global viewpoint, the phenomena safe jump and chaos both imply working unreliability of engineering systems, as they are sensitive to the initial state. An archetypal SD oscillator was proposed in [23] whose abundant dynamics such as bistability and chaos was exhibited. On this basis, Yang and Zhou [24] proposed several feedback-control approaches to control its chaotic vibrations and numerically evaluated the effect of the control. Han and his coworkers proposed another nonlinear oscillator with strong irrational nonlinearities, i.e., the coupled SD oscillator, employed Melnikov method to obtain the analytical criteria of chaotic thresholds, and presented chaos via numerical simulation [25]. Based on SD oscillator, Yang et al. [26] designed a multidirectional multistable device (MMD) including displacement sensors, hydraulic shakers, and oscilloscopes to observe multistability of the oscillator which may be adopted for energy harvesting from ultra-low-frequency vibration sources. Zhang et al. [27] also designed an SD-type isolation system with quasi-zero stiffness and found that the phenomenon bursting oscillation of the system can be triggered by the transition of the stable motion from equilibrium to the limit-cycle oscillation and a jump of solution branch. However, due to the difficulties in dealing with the geometrically nonlinear terms of these oscillatory systems, the mechanism of their dynamics related to initial-state sensitivity is still unclear yet.In this work, we consider a class of nonlinear oscillator with irrational terms and study the mechanism of its various initial-sensitive dynamical behaviors. The study is organized as follows. In Section2, the dynamical model of a typical linkage-slider structure is briefly introduced, and its static bifurcation is then discussed. In Section 3, the average method enabling to describe the primary resonance of the oscillator is illustrated. The phenomenon safe jump is shown via classification of the basins of attraction of multiple periodic attractors. In Section 4, global bifurcation is analyzed. Further numerical studies address geometric aspects such as occurrence of heteroclinic bifurcation and the induced complex dynamics. Finally, Section 5 contains conclusions.


## 2. Dynamical Model

We consider a typical linkage-slider structure [12, 15] whose simplified diagram is shown in Figure 1. The oscillation of the slider is governed by the following equation:(1)mx¨+cx˙+kx−TxL2−x2=F0cosΩt.in which c represents the damping coefficient, k the coefficient of linear stiffness, T the tensile force loading on the ball B, L the length of the connecting rod between the slider A and the ball B, F0 the amplitude of the external excitation on the slider A, and Ω the frequency of the external excitation. The displacement of the slider A should satisfy x≤L. Once the displacement x reaches ±L, the connecting rod will stay horizontally and then cannot move, meaning an unwanted phenomenon of the structure, namely, locking. System (1) becomes dimensionless by letting u=x/L, ω02=k/m, ξ=c/2mω0, τ=ω0t, ω=Ω/ω0, α=T/kL, and f0=F0/kL, expressed as(2)u¨+2ξu˙+u−αu1−u2=f0cosωτ.Figure 1
Simplified oscillator of a linkage-slider structure.Supposing that the tensile forceT in Figure 1 is positive, in the following part, we discuss the case for α>0.The unperturbed system of the dimensionless system (2) can be expressed as(3)u˙=v,v˙=−u+αu1−u2.Accordingly, the number and stability of equilibria of system (2) can be changed when α increases. If 0<α<1, there will be three equilibria, namely, 0,0, 1−α2,0, and −1−α2,0; the origin 0,0 is stable, while the other two equilibria ±1−α2,0 are unstable. If α>1, the only one equilibrium 0,0 will be unstable, which shows that there will be an unbounded escape from the origin; thus, the structure undergoes locking whatever the initial conditions are, which is so-called static locking. The pitchfork bifurcation along α is shown in Figure 2. In order to address more abundant dynamical behaviors of system (2), in what follows, we suppose 0<α<1.Figure 2
Static bifurcation of system (2) with α.


## 3. Multistability and Frequency Jump

Since oscillators are well known to have the complex dynamical behavior in the neighborhood of the main resonance, in this section we follow our previous work [12] and consider that the motion appears to be in the vicinity of the primary resonance. The approximate periodic solution is supposed as u=acosωt−θ. By using the average method, we obtain the following equations of the amplitude and phase(4)a˙=−12πω0∫02π2ω0ξaωsinϕ+δω02acosϕ+αacosϕ1−a2cos2ϕ+f0cosϕ+θsinϕdϕ,θ˙=−12πaω0∫02π2ω0ξaωsinϕ+δω02acosϕ+αacosϕ1−a2cos2ϕ+f0cosϕ+θcosϕdϕ.Letting the right end of the equation above be zero, we can solve the frequencyω and amplitude a in the equation below(5)2ξωa−f0sinθ=0,πω2−1a2+πf0acosθ−4αEa2+4αKa2=0.where the functions Ka2 and Ea2 are the complete elliptic integrals of the first type and the second type, respectively. Conveniently, the numerical results of these two functions can be directly computed by using the symbolic software. Eliminating the parameter θ in equation (5), we have(6)1πa24αEa2−4αKa2−πω2−1a22+4ξ2a2ω2=f02.The theoretical amplitudea can be obtained by solving equation (6) semi-numerically. According to the Jacobian matrix of the linearized equation of the average equation, the corresponding characteristic equation of the periodic solution is(7)λ2+12α3−2ξλ+η=0,where η=−α3ξ+αKa2−Ea24αKa2−4αEa2+πδωa2/π2ω2α3. Obviously, the stability of the periodic solutions can be ascertained by the above equation. It follows that the periodic solution expressed by u=acosωt−θ will be asymptotically stable if η>0.The response curve expressed by equation (6) for ξ=0.001 and α=0.1 is shown in Figure 3. Under different values of the excitation amplitude f0, the variation of the amplitude and stability of periodic solutions with the excitation frequency ω is presented. When f0 is lower, for instance, f0=0.1 or f0=0.3 (see the red or green curve, respectively), there is only one periodic attractor whose amplitude changes continuously with the increase of ω. Comparatively, under a higher f0 like f0=0.5, there is a large range of the excitation frequency that induces locking (see the cyan dashing line for a=1 in Figure 3). Via MATLAB, we apply the 4th Runge–Kutta approach to simulate the numerical solutions of the system so as to verify the validity of our theoretical results. Here, we emphasize that in a certain range of the excitation frequency, different initial conditions of system (2) may lead to different attractors. For example, under the same values of system parameters, i.e., ω=1 and f0=0.1, two different initial states u0,u˙0=0,0 and u0,u˙0=−0.1,0.9 lead periodic oscillation and locking, respectively, as shown in Figure 4(a). This type of locking is termed dynamical locking as the structure does not necessarily undergo locking. In other words, whether locking will occur depends on the initial conditions. Comparatively, when the excitation amplitude is higher, namely, f0=0.5, the structure will undergo static locking (Figure 4(b)).Figure 3
Variation of the amplitude of periodic solutions with the frequency of the excitation whenα=0.1.Figure 4
Time histories of system (2) under different values of f0 when α=0.1 and ω=1. (a) f0=0.1. (b) f0=0.5.
(a)(b)As under certain values of system parameters dynamical locking coexisting with the periodic attractor (Figure4(a)), it is essential for us to classify their basins of attraction. As we know, locking is an unwanted phenomenon for the structure. Its basin of attraction can be understood as a dangerous basin. In contrast, the other initial conditions will lead to bounded solutions; thus, their union is defined as safe basin [28]. Locking instability of the structure can be indicated by fractal erosion of safe basin of system (2), as a small perturbation of initial state for bounded motions may cause locking. We apply the generalized point mapping method to system (2). The domain D=u0≤1,u˙0<2 is partitioned into 200×400 array points. And the trajectories are generated to determine the image of safe basin by fourth-order Runge–Kutta algorithm with one-step mapping length T=0.01. We suppose that if a trajectory satisfies uT<1 within 106 excited circles, then its initial state will be considered safe, thus marked in black. The basin of locking will be marked by white dots. Under the same values of the system parameters ξ, ω, and α of Figure 4(a), one can easily observe the evolution of safe basins with the increase of the excitation amplitude f0 in Figure 5 where the blue “+” indicates the origin 0,0. When f0 increases, the area of safe basin will become smaller, and the fractal of basin boundary will be more obvious. Specifically, when f0=0.42 (Figure 5(d)), the structure may be locked even if the initial conditions are chosen in the neighborhood of the origin. To be worse, when f0=0.5, static locking will occur, for the whole initial domain becomes white (Figure 5(f)). It also agrees with the analytical results in Figure 3. It follows from Figure 5 that the increase of the excitation amplitude f0 will lead to locking instability or even static locking.Figure 5
Evolution of safe basins of system (2) with the increase of f0 when ω=1 and α=0.1. (a) f0=0.1. (b) f0=0.3. (c) f0=0.35. (d) f0=0.42. (e) f0=0.48. (f)f0=0.5.
(a)(b)(c)(d)(e)(f)Given thatf0=0.1, the response curves of equation (5) for different values of the parameter α are shown in Figure 6 where the stable and unstable branches are expressed by solid and dashed curves, respectively. At α=0.1, as shown in the red curve of Figure 6, the unique periodic solution keeps stable, and its amplitude varies continuously with the excitation frequency ω. However, under a higher α, the dynamics will be totally different and sensitive to ω. For instance, when α=0.3, there are two critical values of ω for Hopf bifurcation of system (2), i.e., 0.69 and 0.72, implying that when the excitation frequency ω ranges from 0.69 to 0.72, there will be safe jump between bistable periodic attractors; if ω is beyond the small interval [0.69, 0.72], there will be only one periodic solution which keeps stable and changes continuously with ω, the similar as the case of α=0.1. This evolution of dynamical behaviors with ω can also be verified by the time history diagrams of Figures 7(a) and 7(b). In Figure 7(a), the only periodic response coexists with dynamical locking under the same values of parameters and different initial conditions. In Figure 7(b), when ω=0.7, i.e., within the interval [0.69, 0.72], three different initial conditions lead to bistable periodic responses and static locking, respectively. Under a higher value of α, i.e., α=0.7 (see the purple curves in Figure 6), there are two joints between the dashing curves and solid ones where ω is 0.3 and 0.35, respectively. When α = 0.7, if ω is within the interval [0.30,0.35], static locking will occur (see the horizontal dashing line a¯=1 in Figure 6 as well as the time history diagram in Figure 7(c)), implying that the structure will be definitely in locking despite of the change of initial conditions; if ω is beyond the interval, dynamical locking and periodic response coexist (Figure 7(d)).Figure 6
Amplitude-frequency responses under positive values ofα for f0=0.1.Figure 7
Time histories of system (2) under different values of α and ω. (a) α=0.3, ω=0.35. (b) α=0.3, ω=0.7. (c) α=0.7, ω=0.35. (d) α=0.7, ω=0.7.
(a)(b)(c)(d)The coexistence of periodic attractors and dynamical locking, as shown in Figures7(b) and 7(d), may lead to the loss of global integrity. Hence, we illustrate the evolution of basins of attraction of system (2) with the excitation amplitude f0 in Figure 8 where the black, red, and white regions represent the basins of attraction of the low-amplitude periodic attractor, the high-amplitude one, and locking, respectively. In fact, in Figures 8(a), 8(b), and 8(c) where the two periodic attractors coexist, with the increase of f0, the red region will expand, while the black region will shrink with a clear boundary, implying that safe jump may occur when the initial conditions are chosen near the boundary of basin of attraction of the low-amplitude periodic attractor. In Figure 8(c), the low-amplitude attractor has very all basin of attraction and becomes so-called rare attractor [29]. When f0 continues to increase, the basin of attraction of the low-amplitude periodic attractor will disappear, and the basin boundary of the high-amplitude periodic attractor will be more fractal (Figures 8(d) and 8(e)), showing locking instability. For example, when f0=0.55 (Figure 8(e)), the structure will undergo locking even if the initial conditions are chosen in the neighborhood of the origin. Once f0 increases to 0.65 (Figure 8(f)), the whole initial-condition plane becomes white, which means static locking of the structure.Figure 8
Evolution of basins of attraction of system (2) with the increase of f0 when α=0.3 and ω=0.7. (a) f0=0.1. (b) f0=0.28. (c) f0=0.38. (d) f0=0.45. (e) f0=0.55. (f) f0=0.65.
(a)(b)(c)(d)(e)(f)In this section, it is found that when0<α<1, there are initial-sensitive dynamical behaviors of system (2) such as safe jump and locking instability. The former can be attributed to local bifurcation of periodic solutions and the initial conditions of system (2). We will further discuss the latter in the next section.


## 4. Complex Dynamics Induced by Global Bifurcation

### 4.1. Heteroclinic Orbits of the Unperturbed System

The unperturbed system is given in equation (3) whose Hamiltonian is(8)Hu,v=12v2+12u2+12α2+α1−u2−12.Its corresponding potential energy function is(9)Vu=12u2+α1−u2−12+12α2.According to equations (8) and (9), the existence of potential wells, the shape, the location of possible potential wells, and the number of equilibria of the system depend on the parameter α. If 0<α<1, there will be three equilibria among which the origin is a center, and the other two equilibria P1−α2,0 and Q−1−α2,0 are symmetrical saddle points on the two heteroclinic orbits containing a single potential well and satisfying that(10)Hu,v=12v2+12u2+α1−u2=α.Givenα=0.5, the potential energy and unperturbed orbits are shown in Figures 9(a) and 9(b), respectively. Obviously, within the region surrounded by the heteroclinic orbits (see the bold curves in Figure 9(b)), the trajectories around the origin are closed and elliptic.Figure 9
Diagrams for the unperturbed system (15) when α=0.5 and ω0=1. (a) Potential energy diagram. (b) Phase diagram of the unperturbed system.
(a)(b)To analyze the necessary conditions for global bifurcation of a nonautonomous system, the most common way is to apply the Melnikov method. However, considering the irrational nonlinear terms in equation (10), the heteroclinic orbits of system (7) cannot be expressed by the time τ explicitly. Hence, the Melnikov method cannot be applied directly. We have to propose a new variable ϕ to express these two heteroclinic orbits and time τ explicitly so as to apply the Melnikov method then [28]. Here, the variable ϕ is assumed to satisfy that(11)dϕdτ=ΦϕΦϕ+2π=Φϕ.Additionally, with the increase ofϕ, the heteroclinic orbits go from one saddle point to the other and complete a cycle from 0 to 2π, i.e.,(12)ϕ−∞=0,ϕ+∞=π.Based on the unperturbed system (3), we can further construct heteroclinic orbits ±uhϕ as [30](13)uhϕ=1−α2cosϕ.and have(14)12Φuh′2=−∫uh0uhϕu−αu1−u2du.Substituting equations (13) into (14), we obtain that(15)Φ=±1+3α2−cos2ϕ+22α1+α2−cos2ϕ+α2cos2ϕ1−α21−cos2ϕ.According to equations (13) and (15), the two heteroclinic orbits can be expressed analytically as(16)uhϕ=1−α2cosϕ,vhϕ=∓1+3α2−cos2ϕ+22α1+α2−1−α2cos2ϕ2.Substituting equations (15) into (11), we get(17)τ=±2α−1+αAϕG1φ,4α/1+αAϕ+1+αAϕG2φ,4α/1+αAϕα1+αAϕ,where φ=1/2arcsin1−α2cosϕ, Aϕ=2α+1−αsin2ϕ, and the functions G1 and G2 are the incomplete elliptic integrals of the first type and the second type, respectively. The heteroclinic orbits expressed by equation (10) and by the explicit equations (16) are in excellent agreement (see the bold curves and the red “×” in Figure 9(b)), illustrating that the new expression of the orbits can be used to realize the identity transformation.


### 4.2. Necessary Conditions for Heteroclinic Bifurcation

Now we can use the Melnikov method to discuss necessary conditions for heteroclinic bifurcation. Substituting equations (16) and (17) into the Melnikov function of system (2), we have(18)M±τ0=−43ξ1−α2I1+2f01−α2cosI2−ωτ0,where(19)I1=2+1−ααG1π4+arcsinα2,4α1+α2−1+α2α1−αG2π4+arcsinα2,4α1+α2,I2=ω1+3α22BαG112arcsin1−α22,8αB−ωB2αG212arcsin1−α22,8αB,B=1+α1+3α.For the nonautonomous system (2), if there is a τ0 satisfying(20)M±τ0=0,M˙±τ0≠0.then the Melnikov function (18) will have a simple equilibrium. It is a necessary condition for heteroclinic bifurcation of system (2). Obviously, equation (18) may have simple equilibriums only if(21)2f01−α2>43ξ1−α2I1,namely,(22)f0>f0cri=23ξ1−α2I1,where f0cri is the threshold for heteroclinic bifurcation of system (2). The variation of the threshold f0cri with the parameter α under ω=1 is shown in Figure 10(a), while the change of f0cri with the frequency ω under a fixed value of α is shown in Figure 10(b). According to Figure 10(a), the change of f0cri with α is nonmonotonical. The threshold f0cri decreases with the increase of α when α≤0.23, but increases with α when α>0.23. To be different, the threshold f0cri increases monotonically with the frequency ω (Figure 10(b)).Figure 10
Variation of the thresholds off0 for heteroclinic bifurcation with system parameters (a) ω=1. (b) α=0.1.
(a)(b)


### 4.3. Numerical Examples

Global bifurcation often leads to the erosion of safe basin [28] and chaos [12, 18]. We wonder if the heteroclinic bifurcation in system (2) may cause the two types of initial-sensitive phenomena simultaneously. Thus, in this subsection, we discuss these typical complex dynamical behaviors numerically. The numerical approaches and the setting of the initial conditions are the same as in Section 3.First, we obtain numerically minimum values off0 for locking instability under different values of the system parameters ω and α (see the red “×” in Figure 10). As mentioned in Section 3, locking instability can be depicted by fractal erosion of safe basin of system (2).It follows from the comparison of theoretical and numerical results in Figure10 that the numerical results for critical values of f0 that induces locking instability totally agree with the analytical ones, showing that the heteroclinic bifurcation of system (2) can lead to locking instability. It also reflects the accuracy of our analysis.Second, the evolution of safe basin of system (2) with the variation of the system parameters is shown in Figures 11 and 12 where the blue “+” indicates the position of the origin. According to the theoretical results and numerical ones in Figures 3 and 6, it is clear that when ω=1, 0.1≤α≤0.28, and 0.24≤f0≤0.5, system (2) will undergo a unique periodic motion or locking. It means that the safe basin in Figure 11 is the basin of attraction of the only periodic attractor of system (2). In Figures 11 and 12, as the parameter f0 increases, the safe basin will become fractal, and its area will be reduced, implying that locking instability becomes more and more obvious. For example, in Figures 11(c) and 11(g), the basin boundary is obviously fractal. When f0 continues to increase (Figures 11(d), 11(h), 12(d), and 12(h)), the whole initial-condition plane is white, verifying static locking predicted in Figure 3. According to Figure 11, the area of safe basin firstly expands and then reduces with the increase of α, which also agrees with the analytical trend of f0cri in Figure 10(a). Besides, as shown in each column of Figure 12, the increase of ω can reduce the basin erosion, i.e., locking instability, which also verifies the predicted trend of f0cri in Figure 10(b).Figure 11
Evolution of safe basins of system (2) with the variation of α and f0 when ω=1. (a) α=0.1, f0=0.24. (b) α=0.1, f0=0.26. (c) α=0.1, f0=0.4. (d) α=0.1, f0=0.5. (e) α=0.2, f0=0.24. (f) α=0.2, f0=0.26. (g) α=0.2, f0=0.4. (h) α=0.2, f0=0.5. (i) α=0.28, f0=0.24. (j) α=0.28, f0=0.26. (k) α=0.28, f0=0.4. (l) α=0.28, f0=0.5.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Figure 12
Evolution of safe basins of system (2) with the variation of ω and f0 when α=0.1. (a) ω=0.9, f0=0.21. (b) ω=0.9, f0=0.23. (c) ω=0.9, f0=0.35. (d) ω=0.9, f0=0.75. (e) ω=1.1, f0=0.21. (f) ω=1.1, f0=0.23. (g) ω=1.1, f0=0.35. (h) ω=1.1, f0=0.75. (i) ω=1.2, f0=0.21. (j) ω=1.2, f0=0.23. (k) ω=1.2, f0=0.35. (l) ω=1.2, f0=0.75.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Finally, we observe if the chaotic motion will occur in system (2) when the heteroclinic bifurcation of system (2) is induced. According to equation (22), we can calculate that the critical value of f0 for heteroclinic bifurcation is 0.69. Now fixing the values of the parameters α and ω, i.e., α=0.1 and ω=1.2, and varying the values of the parameter f0 from 0 to 1, the bifurcation diagram of system (2) in Poincare map is shown in Figure 13. When f0 is less than 0.69 (see the left side of the vertical dashing line in Figure 13), system (2) will be in periodic motion. In case that f0 exceeds 0.69, one can easily observe from Figure 13 that system (2) evolves chaotic motion through the routine of the period-doubling bifurcation, illustrating that heteroclinic bifurcation of system (2) also leads to chaos.Figure 13
Bifurcation diagram of system (2) with the increase of f0 when α=0.1 and ω=1.2.


## 4.1. Heteroclinic Orbits of the Unperturbed System

The unperturbed system is given in equation (3) whose Hamiltonian is(8)Hu,v=12v2+12u2+12α2+α1−u2−12.Its corresponding potential energy function is(9)Vu=12u2+α1−u2−12+12α2.According to equations (8) and (9), the existence of potential wells, the shape, the location of possible potential wells, and the number of equilibria of the system depend on the parameter α. If 0<α<1, there will be three equilibria among which the origin is a center, and the other two equilibria P1−α2,0 and Q−1−α2,0 are symmetrical saddle points on the two heteroclinic orbits containing a single potential well and satisfying that(10)Hu,v=12v2+12u2+α1−u2=α.Givenα=0.5, the potential energy and unperturbed orbits are shown in Figures 9(a) and 9(b), respectively. Obviously, within the region surrounded by the heteroclinic orbits (see the bold curves in Figure 9(b)), the trajectories around the origin are closed and elliptic.Figure 9
Diagrams for the unperturbed system (15) when α=0.5 and ω0=1. (a) Potential energy diagram. (b) Phase diagram of the unperturbed system.
(a)(b)To analyze the necessary conditions for global bifurcation of a nonautonomous system, the most common way is to apply the Melnikov method. However, considering the irrational nonlinear terms in equation (10), the heteroclinic orbits of system (7) cannot be expressed by the time τ explicitly. Hence, the Melnikov method cannot be applied directly. We have to propose a new variable ϕ to express these two heteroclinic orbits and time τ explicitly so as to apply the Melnikov method then [28]. Here, the variable ϕ is assumed to satisfy that(11)dϕdτ=ΦϕΦϕ+2π=Φϕ.Additionally, with the increase ofϕ, the heteroclinic orbits go from one saddle point to the other and complete a cycle from 0 to 2π, i.e.,(12)ϕ−∞=0,ϕ+∞=π.Based on the unperturbed system (3), we can further construct heteroclinic orbits ±uhϕ as [30](13)uhϕ=1−α2cosϕ.and have(14)12Φuh′2=−∫uh0uhϕu−αu1−u2du.Substituting equations (13) into (14), we obtain that(15)Φ=±1+3α2−cos2ϕ+22α1+α2−cos2ϕ+α2cos2ϕ1−α21−cos2ϕ.According to equations (13) and (15), the two heteroclinic orbits can be expressed analytically as(16)uhϕ=1−α2cosϕ,vhϕ=∓1+3α2−cos2ϕ+22α1+α2−1−α2cos2ϕ2.Substituting equations (15) into (11), we get(17)τ=±2α−1+αAϕG1φ,4α/1+αAϕ+1+αAϕG2φ,4α/1+αAϕα1+αAϕ,where φ=1/2arcsin1−α2cosϕ, Aϕ=2α+1−αsin2ϕ, and the functions G1 and G2 are the incomplete elliptic integrals of the first type and the second type, respectively. The heteroclinic orbits expressed by equation (10) and by the explicit equations (16) are in excellent agreement (see the bold curves and the red “×” in Figure 9(b)), illustrating that the new expression of the orbits can be used to realize the identity transformation.


## 4.2. Necessary Conditions for Heteroclinic Bifurcation

Now we can use the Melnikov method to discuss necessary conditions for heteroclinic bifurcation. Substituting equations (16) and (17) into the Melnikov function of system (2), we have(18)M±τ0=−43ξ1−α2I1+2f01−α2cosI2−ωτ0,where(19)I1=2+1−ααG1π4+arcsinα2,4α1+α2−1+α2α1−αG2π4+arcsinα2,4α1+α2,I2=ω1+3α22BαG112arcsin1−α22,8αB−ωB2αG212arcsin1−α22,8αB,B=1+α1+3α.For the nonautonomous system (2), if there is a τ0 satisfying(20)M±τ0=0,M˙±τ0≠0.then the Melnikov function (18) will have a simple equilibrium. It is a necessary condition for heteroclinic bifurcation of system (2). Obviously, equation (18) may have simple equilibriums only if(21)2f01−α2>43ξ1−α2I1,namely,(22)f0>f0cri=23ξ1−α2I1,where f0cri is the threshold for heteroclinic bifurcation of system (2). The variation of the threshold f0cri with the parameter α under ω=1 is shown in Figure 10(a), while the change of f0cri with the frequency ω under a fixed value of α is shown in Figure 10(b). According to Figure 10(a), the change of f0cri with α is nonmonotonical. The threshold f0cri decreases with the increase of α when α≤0.23, but increases with α when α>0.23. To be different, the threshold f0cri increases monotonically with the frequency ω (Figure 10(b)).Figure 10
Variation of the thresholds off0 for heteroclinic bifurcation with system parameters (a) ω=1. (b) α=0.1.
(a)(b)


## 4.3. Numerical Examples

Global bifurcation often leads to the erosion of safe basin [28] and chaos [12, 18]. We wonder if the heteroclinic bifurcation in system (2) may cause the two types of initial-sensitive phenomena simultaneously. Thus, in this subsection, we discuss these typical complex dynamical behaviors numerically. The numerical approaches and the setting of the initial conditions are the same as in Section 3.First, we obtain numerically minimum values off0 for locking instability under different values of the system parameters ω and α (see the red “×” in Figure 10). As mentioned in Section 3, locking instability can be depicted by fractal erosion of safe basin of system (2).It follows from the comparison of theoretical and numerical results in Figure10 that the numerical results for critical values of f0 that induces locking instability totally agree with the analytical ones, showing that the heteroclinic bifurcation of system (2) can lead to locking instability. It also reflects the accuracy of our analysis.Second, the evolution of safe basin of system (2) with the variation of the system parameters is shown in Figures 11 and 12 where the blue “+” indicates the position of the origin. According to the theoretical results and numerical ones in Figures 3 and 6, it is clear that when ω=1, 0.1≤α≤0.28, and 0.24≤f0≤0.5, system (2) will undergo a unique periodic motion or locking. It means that the safe basin in Figure 11 is the basin of attraction of the only periodic attractor of system (2). In Figures 11 and 12, as the parameter f0 increases, the safe basin will become fractal, and its area will be reduced, implying that locking instability becomes more and more obvious. For example, in Figures 11(c) and 11(g), the basin boundary is obviously fractal. When f0 continues to increase (Figures 11(d), 11(h), 12(d), and 12(h)), the whole initial-condition plane is white, verifying static locking predicted in Figure 3. According to Figure 11, the area of safe basin firstly expands and then reduces with the increase of α, which also agrees with the analytical trend of f0cri in Figure 10(a). Besides, as shown in each column of Figure 12, the increase of ω can reduce the basin erosion, i.e., locking instability, which also verifies the predicted trend of f0cri in Figure 10(b).Figure 11
Evolution of safe basins of system (2) with the variation of α and f0 when ω=1. (a) α=0.1, f0=0.24. (b) α=0.1, f0=0.26. (c) α=0.1, f0=0.4. (d) α=0.1, f0=0.5. (e) α=0.2, f0=0.24. (f) α=0.2, f0=0.26. (g) α=0.2, f0=0.4. (h) α=0.2, f0=0.5. (i) α=0.28, f0=0.24. (j) α=0.28, f0=0.26. (k) α=0.28, f0=0.4. (l) α=0.28, f0=0.5.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Figure 12
Evolution of safe basins of system (2) with the variation of ω and f0 when α=0.1. (a) ω=0.9, f0=0.21. (b) ω=0.9, f0=0.23. (c) ω=0.9, f0=0.35. (d) ω=0.9, f0=0.75. (e) ω=1.1, f0=0.21. (f) ω=1.1, f0=0.23. (g) ω=1.1, f0=0.35. (h) ω=1.1, f0=0.75. (i) ω=1.2, f0=0.21. (j) ω=1.2, f0=0.23. (k) ω=1.2, f0=0.35. (l) ω=1.2, f0=0.75.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Finally, we observe if the chaotic motion will occur in system (2) when the heteroclinic bifurcation of system (2) is induced. According to equation (22), we can calculate that the critical value of f0 for heteroclinic bifurcation is 0.69. Now fixing the values of the parameters α and ω, i.e., α=0.1 and ω=1.2, and varying the values of the parameter f0 from 0 to 1, the bifurcation diagram of system (2) in Poincare map is shown in Figure 13. When f0 is less than 0.69 (see the left side of the vertical dashing line in Figure 13), system (2) will be in periodic motion. In case that f0 exceeds 0.69, one can easily observe from Figure 13 that system (2) evolves chaotic motion through the routine of the period-doubling bifurcation, illustrating that heteroclinic bifurcation of system (2) also leads to chaos.Figure 13
Bifurcation diagram of system (2) with the increase of f0 when α=0.1 and ω=1.2.


## 5. Conclusions

In order to investigate the mechanism of initial-sensitive dynamical behaviors of an oscillatory system with geometrical nonlinearities, we consider a typical linkage-slider structure as the dynamical model. Its complex dynamics such as safe jump and locking instability, and chaos is discussed in detail. The study shows that the dynamical system behaves as pitchfork bifurcation, Hopf bifurcation, and heteroclinic bifurcation as the system parameters such as tensile force, the excitation amplitude, and the excitation frequency vary. The agreement of analytical results and numerical ones verifies the accuracy of our theoretical prediction. The main conclusions are presented as follows.(1)
Static locking of the structure may occur under a heavy tensile force, which is attributed to pitchfork bifurcation(2)
When the tensile force and the excitation amplitude are low, within a certain range of the excitation frequency, safe jump between bistable periodic attractors can be induced by Hopf bifurcation and a slight perturbation of initial conditions(3)
Due to heteroclinic bifurcation of the oscillatory system, locking instability and chaos of the structure may occur. The increase of the excitation amplitude will induce locking instability, chaos, or even static locking, while the increase of the excitation frequency can reduce locking instability. The variation of the tensile force can lead to locking instability, but the extent of locking instability does not change monotonically with it.This study provides some references in the design and applications of mechanical oscillators with the irrational nonlinearities. Nevertheless, the case for the negative tensile force has not been discussed yet where different bifurcations and rich dynamics may be induced. Further, referring to [26], we may construct experiments to observe the initial-sensitive phenomena such as safe jump, locking instability, and chaos, which will be included in our future work.



---

*Source: 6472678-2022-06-30.xml*