On the primary resonance of Duffing oscillator: a novel approximate algebraic equation for solving jump points and a robust criterion for global stability

Duffing oscillator is a classical nonlinear system widely used in various fields of science and engineering. In this paper, jump points and global stability of primary resonance of Duffing system, including softening and hardening, are further investigated. A novel higher order polynomial equation in terms of amplitudes is determined to calculate the jump points of primary frequency response curve (FRC) by using a first-order approximate HBM. The jump points determined by polynomial equation are in good with the results of direct numerical simulation, and those given in other previous studies under the same system parameters. Further, explicit analytical expressions of the boundary curves (BCs), that can be used to quickly determine the unstable response region of the primary FRC, is deduced by using Floquet theory. The BCs can also be considered as a set of jumping points of the FRCs under different excitations. According to a necessary condition of separation of FRC and BCs, the global stability criterion of FRC is transformed into the relationship between BCs and the extreme points of FRC. Noteworthily, the global stability criteria of FRC for both softening and hardening systems, that are analytical expressions composed of damping ratio, nonlinear coefficient and external excitation, can be determined. The criterion is provided a novel quantitative analysis method for the design of a Duffing nonlinear system to avoid or utilize the jump characteristic.


Introduction
Recent years, Duffing oscillator, an old-aged and classical nonlinear system, has been widely used to describe approximate models of the system in engineering, such as quasi-zero stiffness vibration isolation [1][2][3][4], vibration energy harvesting [5][6][7]22], and micro-and nanoelectromechanical systems [8,9]. Duffing oscillator exhibits an enormous range of well-known behaviours in nonlinear dynamical systems, including homoclinic bifurcations, catastrophe, chaos, super-and sub-harmonic, etc [10]. Bifurcation of the frequency response curve (FRC) also is called the jump phenomenon, which is closely related to the stability of periodic response of the system. In this paper, a novel higher order polynomial equation in terms of amplitudes is determined the jump points of primary FRC by using a first-order approximate harmonic balance method (HBM). Furthermore, global stability criteria of the FRC for both softening and hardening systems, which are analytical expressions composed of damping ratio, nonlinear coefficient and external excitation, can be determined by using Floquet theory.
A large number of approximate analytical methods, including classical perturbation methods (PMs) [11], multi-scale method (MSM) [12], HBM [13] and improved harmonic balance methods [14,15], homotopy method [16] and improved homotopy methods [17,18], as well as numerical methods, including time simulation, periodic extension and shooting method, have been applied to the solution of Duffing equation and analysis of its dynamic behaviours. Nayfeh and Mook [12] used MSM to determine the jump-down frequency and peak amplitude of the response. Kevorkian and Cole [19] determined the jump up and down frequencies and its corresponding amplitudes for a hardening system by the same method. Magnus [20] used HBM to determine the frequency of the jump-down point and the peak amplitude of a hardening system but did not consider the jump-up point. Brennan and Kovacic [21] determine the jump-down frequency and the corresponding amplitudes by HBM. They also summarised the approximate expressions of the jump-up point and jump thresholds given by PMs. Mallik [10] pointed out that FRC has vertical slope at the jump point, where dω/dZ=0; here, ω and Z represent the frequency and amplitude, respectively. Unfortunately, they did not elaborate the solution of the jump points. Recently, some literatures used these expressions [21] in the engineering application researches, such as nonlinear energy-harvesting devices [6,22], Duffing-type suspension system [23], high-static-low-dynamic-stiffness isolator [24]. Ho and Lang [25] determined the output frequency response function to examine analytically the jump points of a Duffing oscillator. In the approximate analytical method [12,[19][20][21][22][23][24][25], it is generally assumed that the damping of the system is small, so that the damping term above the second order in frequency response function (FRF) can be ignored to obtain analytical expressions. Therefore, the jump-down point is approximated to the maximum response of FRF, and the jumpup point given by PMs is considered to be a damping independent expression. However, with the increase of damping, we find that the higher-order damping term in FRF still affects the value of the jump point.
In numerical and hybrid analytical-numerical methods, Peleg [26], studying a hardening isolation system with friction and viscous damping, proposed an implicit amplitude-frequency equation. The jump frequencies and the corresponding amplitudes are determined by locating the vertical tangents to the FRC. Friswell and Penny [27] used a numerical approach based on Newton's method, including terms up to the ninth harmonic, to compute the "exact" frequencies of the jump points. Worden [28] used a first-order HBM expansion and set the discriminant of a cubic polynomial in the square of amplitude equal to zero to obtain a higher-order polynomial equation in terms of jump frequency. Compared with the "exact" values of Friswell [27], the results show that the numerical solutions of the higher order polynomial proposed by Worden [28] is accurate enough. It is also noteworthy that Malatkar and Nayfeh [29] solve the jump points by constructing the Sylvester resultant of the first and second derivatives of FRF. A sixth order polynomial equation on jump frequencies determined by Sylvester resultant can provide accurate jump values, but it needs to be solved by a simple numerical method. They also given the expressions of jump frequency and jump threshold when the damping coefficient is equal to zero, which are approximately applied to the weak damping system. Numerical methods are accurate for solving values of jump points, but the disadvantage is that the solutions may converge slowly or to incorrect values.
Inspired by the Worden [28] and Malatkar [29], and using the property that FRC has a vertical slope at the jump points, a quintic polynomial equation in term of the square of jump amplitudes is derived in this paper, which is more concise than the equations in references [28,29]. The corresponding jump frequencies can be obtained by substituting the jump amplitudes into FRC. In order to make the polynomial equation converge to an accurate value quickly, we use the approximate analytical formulas summarized in reference [21] to calculate the initial values of numerical simulation. The jump points determined by the polynomial equation given in this paper are in good agreement with the results of direct numerical simulation, and those given by Friswell [27] and Worden [28] under the same system parameters.
In addition, we use the Floquet theory to further study the global stability of the primary FRC of the system to determine the conditions for jump occurrence. Conversely, whether the FRC jumps at a certain point is determined by the stability of the response at that point. Lyapunov [30] and Floquet [31] theories are the commonly used methods to study the response stability of nonlinear systems, which mainly use the roots of characteristic equation to determine the stability of various systems [32][33][34][35].
In this study, the characteristic equation of the system is further expanded to obtain the explicit analytical expression of BCs, which can be used to quickly determine the unstable response region of FRC. The intersections of BCs and primary FRC are the jump points. More importantly, we further discuss how to determine the global stability of FRC, that is, the condition that there is no jump point in FRC. We find that the global stability of the system is dominated by the topological relationship between FRC and BCs, and can be transformed into the topological relationship between FRC and the extreme point of BCs. Finally, the global stability criteria in softening and hardening systems are determined, which are analytical expressions composed of damping ratio, nonlinear coefficient and external excitation. The global stability criterion is provided a fast and quantitative analysis method for the design of a Duffing nonlinear systems to avoid or utilize the jump characteristic.

Amplitude-frequency response
The differential equation of Duffing oscillator excited by a harmonic force Acos(ωt+φ) is given by ( ) 3 13 cos mz cz k z k z A t  where m is the equivalent mass of the system, z is the displacement, c is the damping, and k1 and k3 are the linear and cubic stiffness terms, respectively. The non-dimensional form of Eq. (1) can be written as Parameter a denote the equivalent excitation of the system, and γ is called nonlinear coefficient that is the ratio of cubic stiffness to linear stiffness.
A complete FRC includes amplitude-frequency response (AFR) and phase-frequency response (PFR). The solution approach of jump points of AFR is explained here firstly. Applying HBM, assuming a first-order approximate solution in the form of ( ) cos zZ  = , and neglecting the term of cos3Ωτ, the implicit AFR is given as It can be expanded as a quadratic equation for Ω 2 , Then, the explicit AFR derived from Eq. (4) By setting Eq. (5a) and Eq. (5b) to be an equality, a quadratic equation in the form of Z 2 can be determined ( ) Solving Eq. (6) and deleting the negative root, the maximum amplitude of AFR is given as Substituting Eq. (7) into Eq. (5a) or (5b) gives the corresponding frequency as ( ) When damping ratio ζ is small enough, the term ζ 2 can be neglected from Eqs. (7) and (8) to get the same results as those in reference [21]. Eqs. (7) and (8), which are valid for both hardening and softening systems, yield positive real values that satisfy the following inequalities: and also ( ) For hardening system with periodic response (γ>0, 0≤ζ≤1), Eqs. (9a), (9b) and (9c) are always satisfied. However, it can be deduced from Eqs. (9a), (9b) and (9c) that the external excitation of softening system with periodic response (γ<0, 0≤ζ≤1), should be less than a maximum value, which is given as If a is greater than amax, Ωm and Zm will become complex, which mean that the potential energy cannot overcome the external excitation, so that the system no longer has periodic response. Substituting Eq. (6) into Eq. (5a) or Eq. (5b), the damping backbone curve which represents the damping free oscillation of system is given by Fig. 1 shows the AFR of Duffing system with damping ratio ζ=0.1 and excitation a=1, as well as their damping backbone curves, maximum response and jump points. The AFR for hardening system bends to the right, while that for softening system bends to the left. The damping backbone curve intersects with AFR at the maximum response point.

A higher-order polynomial equation for calculating jump points
AFR have vertical slopes at both jump down and up points, that is dZ/dΩ=∞ or dΩ/dZ=0 [10]. Therefore, applying dΩ/dZ=0 to the explicit AFR (Eqs. (5a, b)), an equation with amplitudes of the jump points can be written as Expanding and simplifying Eq. (12), a quintic polynomial equation in terms of the square of jump amplitude is determined as 2 In references [28,29], higher-order polynomial equations on jump frequencies have been given from different techniques, while in this paper, corresponding solutions have been given from the perspective of jump amplitude, and the polynomial equation presented in this paper is more concise. According to Galois theorem, there are no radical solutions for quintic equation and equations with more than five orders. In addition, it is complicated to use elliptic function to solve higher-order equation, which is not involved in this paper. Therefore, Newton's method was used to solve Eq. (13) in this paper. Firstly, Eq. (13) is rewritten as function F(Z). The iterative equation is defined as where i is the number of iterations, Zi+1 is the amplitude after the i-th iteration. The iterative initial values are calculated by the formulas summarized in reference [21], which can make the system converge to accurate values rapidly. The iterative termination condition is set to 1 ii ZZ  + − , where ε is the precision threshold. Then, substituting jump amplitudes into Eq. (5a) or Eq. (5b) yields the corresponding jump frequencies.
The evolution of jump points for hardening system with coefficient γ in the range [γh,3γh] is shown in the Fig.2. Parameter γh called nonlinear threshold for hardening system will be explained in detail in section 3.2. Since different damping and excitation result in different nonlinear coefficient for the system jump, γ is normalized to [0,1] represented by a symbol γn for the convenience of comparison. Symbol Zd denote jump-down amplitude, Ωd denote jump-down amplitude frequency, Zu denote jump-up amplitude and Ωu denote jump-up frequency. Results in Fig.2 show that the values of jump points given by Eq. (13) are in good agreement with the numerical sweep frequency simulation of equation Eq. (2). The jump-up amplitudes are decreasing with increase of γn, while the jump frequencies are increasing. This is due to the increase of the degree of nonlinearity, the AFR of hardening system bends more to the right. When γn and a are constant, jump amplitudes decrease and jump frequencies increase with the damping increases. Similarly, for a softening system, the variation of jump points with γ in the range [1.5γs, γs] is shown in the Fig.3. Parameter γs called nonlinear threshold for softening system will be also defined in detail in section 3.2. Coefficient γ is normalized to [-1,0] represented by the symbol γn. Fig. 3 show that, for the softening system, the jump points given by Eq. (13) also fit well with the direct numerical simulations of Eq. (2). Increasing value γn, the jump frequencies and jump-up amplitudes are increasing, while the jump-down amplitudes are decreasing. When γn and a are set to constant and the damping value increases, the jump frequencies and the corresponding amplitude will decrease. Furthermore, following the system parameters set in references [27,28], where the values m=γ=1, we further verify the accuracy of estimation formulas for the jump points proposed in this paper. The results in Table 1 and 2 show that the values in this paper are completely consistent with those in reference [28], and the percentage error with the "exact" value in reference [27] is basically less than 0.3%.

Calculating for the jump points of PFR
The phase frequency response (PFR) [10] can be written as  Fig. 4. It can be seen from the figure that all the PFRs bend to the right. The larger value γ is, the more obvious the bending is. The phase φm is near π/2 and gradually approaches π/2 with the increase in γ. When γ<γh, there is no jump phenomenon in the PFR, as denoted by the green line. When γ=γh, there is only one jump point, that is, the jump-down point is equal to the jump-up point. When γ>γh, there are two jump points to form an unstable region, which becomes larger as the value of γ increases.

Boundary curve of the steady-state region
Stability of the periodic solutions of the system is calculated by applying Floquet theory as described in [33,34]. To this aim, a small disturbance e(τ) is introduced to the harmonic solutions u(τ) of the equation of motion given in Eq. (2), leading to  (24) Generally, the stability of steady-state motion depends on the eigenvalues of the Jacobian matrix. The eigenvalues can be constructed to a criterion function of stability that is The function Г>0 means that the periodic solution is asymptotically stable, while Г ≤ 0 represents unstable periodic response. Consequently, Г=0 can determine the boundary between the steady and unsteady periodic responses of the system. A quadratic equation is obtained in terms of BCs, given by Eqs. (27a) and (27b), plotted in the frequency-amplitude plane, provides the regions of instability (filled red areas) as shown in Fig. 5 and 6. The region is a function of the damping ratio ζ and nonlinear coefficient γ but it does not depend on the displacement excitation a. When the external excitation level a is low, the AFR curve is low and does not reach the unstable region, which indicates that the AFR is globally stable. While the external excitation a exceeds a certain critical value, part of AFR will cross the BCs to form unstable responses. These intersections between BCs and AFRs are the jump points.  The results in Fig. 7 show that the regions of instability formed by BCs varies with damping ratio ζ. When the nonlinear coefficient γ is constant, the unstable region of AFR decreases with the increase of ζ, whether it is hardening system or softening system. Those unstable regions with high damping ratio are surrounded by the unstable regions with low damping ratio. It should be noted that the jump points of AFRs just falls on its corresponding BCs. The jump-down point of undamped Duffing system is at infinity, so it cannot be plotted in Fig. 7.  When the damping ratio ζ is constant, the unstable region of AFR also decreases with the increase of γ for both hardening system and softening system. Specifically, the amplitudes of the unstable response become smaller with the increase of γ, while the corresponding frequency range has no obvious change.

Global stability criteria of primary AFR
Extreme points of the BCs given by Eqs. (27) can be obtained by equal the discriminant of Eq. (26) to zero, that is ( ) Solving the quadratic Eq. (28) in term of Ω 2 , the frequencies of extreme points are given as where Uh given by Eq. (30a) is used for hardening system and Us given by Eq. (30b) is used for softening system. Considering the condition ζ>0 in this paper, the values of Ωs, Ωh, Us, and Uh, are always positive.
Substituting Us into Eq. (5a) to replace Z, a frequency of AFR for a softening system is determined as Therefore, setting Ωs=Ωfs and simplifying, the excitation threshold of the global stability of a softening system, which also determines whether the system jumps or not, is determined by the following equation  (32) This means that the AFR of a softening system is globally stable when the condition a<as is satisfied. Contrarily, when the condition a≥as is satisfied, the AFR of a softening system will form unstable regions, where the BCs passes through the jump points. Similarly, setting Ωh=Ωfg, the excitation threshold of the global stability of a hardening system, is determined by As shown in Fig. 9, when the nonlinear coefficient γ is constant, the excitation thresholds increase with the increase of damping ratio ζ for both hardening and softening systems. At the same time, when the damping ratio ζ is constant, with the increase of γ, the excitation thresholds will also increase. A hardening system have similar criteria for determining the global stability of its AFR. Therefore, the global stability criteria (jump bifurcation conditions) of the AFRs for softening and hardening systems are summarized in Table 3. Table 3 Global stability criteria of the AFR in Duffing oscillator

Category
Stability of AFR Jump points Criteria