Periodic Solution of the Strongly Nonlinear Asymmetry System with the Dynamic Frequency Method

: In this article, we present a new accurate iterative and asymptotic method to construct analytical periodicsolutions for a strongly nonlinear system, even ifit is not Z 2 -symmetric. This method is applicable not only to a conservative system but also to a non-conservative system with a limit cycle response. Distinct from the general harmonic balance method, it depends on balancing a few trigonometric terms (at most ﬁve terms) in the energy equation of the nonlinear system. According to this iterative approach, the dynamic frequency is a trigonometric function that varies with time t , which represents the inﬂuence of derivatives of the higher harmonic terms in a compact form and leads to a signiﬁcant reduction of calculation workload. Two examples were solved and numerical solutions are presented to illustrate the e ﬀ ectiveness and convenience of the method. Based on the present method, we also outline a modiﬁed energy balance method to further simplify the procedure of higher order computation. Finally, a nonlinear strength index is introduced to automatically identify the strength of nonlinearity and classify the suitable strategies.


Introduction
Nonlinear differential equations describing dynamic behaviors play a major role in mechanics and mathematics. The analysis of the literature devoted to the theory of nonlinear differential equations show that various so-called analytical, functional-analytic, numerical and numerical-analytic methods based upon successive approximations are now extensively studied [1]. The group of numerical methods under the assumption of the existence of periodic solutions gives numerical algorithms for approximate solutions [2]. Compared with the numerical means, the analytical methods are more interesting because such solution allows direct discussion of the effects of related parameters in the qualitative investigation (stability, branching, uniqueness). And the so called numerical-analytic methods [3,4] have been widely used to investigate nonlinear boundary value problems (BVPs).
Among the analytical methods, there are many papers dealing with classical solution methods for nonlinear differential equations which describe the one-degree-of-freedom system. Most of the procedures are of the perturbation type [5,6]: the Krylov-Bogolubov-Mitropolsky (KBM) method [7,8], the averaging solution method [9], the multiple scales method [10], the Lindstedt-Poincare (L-P) perturbation method [11], etc. However, these theories mentioned are based on the presence of a small parameter in the nonlinear governing equation. The requirement of a small parameter assumption greatly restricts the applications of perturbation methods, and it is particularly influential for strongly

The Basic Idea of the Dynamic Frequency Method
Consider the following autonomous nonlinear system ..
where α i,0 , β i,j are the polynomial coefficients; M and K are the integers, M > 1, K ≥ 0 and an overdot represents differentiation with respect to time t.
We assume there exists some equilibrium of positions in Equation (1), around which the system can perform periodic motions. According to the harmonic balance method, the solution to Equation (1) can be expressed in the form a n cos nω 1,0 t, . u = −a 0 ω 1,0 sin ω 1,0 t − N n=2 a n nω 1,0 sin nω 1,0 t.
Based on the expansion formula of trigonometric function, sin nω 1,0 t can be transformed into Γ n (t) sin ω 1,0 t in the case of n ≥ 2, where Γ n (t) is a function of sin ω 1,0 t and cos ω 1,0 t, i.e., Hence Equation (2) can be rewritten as where P n is the coefficient of the nth-order harmonic term.
Since the harmonic balance method is based on the assumption that the fundamental frequency term is dominant in a Fourier representation of the solution, a n is a relatively small quantity compared to a 0 . Thus, by neglecting the higher-order harmonic terms in u, Equation (4) can be represented as where a small nondimensional parameter p has been introduced as a bookkeeping parameter and set equal to unity in the final result, ω 1,0 is the undetermined fundamental frequency and ω 1,i (t) can be assumed as a dynamic frequency which is a function of Γ n (t) in Equation (4). It is worth noting that the high-order harmonic terms are not neglected in the expression of . u for the reason that the amplitude is a slow variable while frequency is a fast variable.
Different from the processure of the harmonic balance method, the unknown variables can be solved in the energy equation of the system. The integral of Equation (1) is . udt, where E * is an average mechanical energy over a whole period. According to Li [41], if the Equation (1) has a stable limit cycle, there should exist E * > 0 such that dE * dt = 0. Substituting Equation (5) into equation (6), performing the integration and collecting the power series of p lower than k for each order, we can obtain the following set of equations • order 2: • order k: It is necessary to balance some trigonometric function terms on both sides of Equation (9) to obtain those unknown variables (E 0 , a 0 , b 0 , ω 1,0 , ω 1,k (t)). Herein, to be distinguished from the general Fourier expansion in HB, the new dynamic frequency approach does not terminate the balancing procedure at any order harmonics cos kω 1,0 t, sin kω 1,0 t but uses ω 1,k (t) as a compact variable to represent all those remaining harmonics to obtain the algebraic equations.
During that course, one needs to balance at most five terms in each order to find those unknown variables for the periodic solutions, no matter what the specific m and n are. That greatly reduces workload and difficulty for quantitative analysis in HB. The basic algorithm is Step 1 : balance the constant term, Step 2 : balance the term of t or sin ω 1,0 t cos ω 1,0 t, Step 3 : balance the term of cos ω 1,0 t, Step 4 : balance the term of sin 2 ω 1,0 t, Step 5 : balance the term of remaining terms.
The derivative relationship between (u, . u) can be permitted through an integration which gives θ(t) as a periodic phase component in the following approximate solution Besides, it is important to note that the integration of ω 1,k (t) in Equation (11) will produce a linear term ∆ 1,k t where the ∆ 1,k can be regarded as the new static frequency component. Then we can introduce the following transformation where ω k is the kth approximate frequency and ∆ 1,k depends on the choice of the order k in ω 1,i (t).
Since the present method is based on the assumption that the fundamental frequency term is dominant in the solution, the assumption must be checked a posteriori, i.e., ∆ 1,k < ∆ 1,k−1 . Replacing the ω 1,0 in all the nonlinear terms of Equation (12) with the frequency ω k , the kth-order approximate solution to the system is determined. To this end, the whole procedure is implemented in an algebraic manipulator, such as MATHEMATICA. In general, the procedure would become increasingly cumbersome as the order goes up. More importantly, the computation results show that the solution up to the second order is fairly accurate, even for strongly nonlinear oscillators. In the following sections, two examples will be presented to demonstrate the effectiveness and applicability of this new method. We also introduce the modified energy balance method to simplify that procedure of higher order computation.

Duffing Oscillators
Duffing equation can be written in the form ..
Based on Equation (7), the first order approximate solution to Equation (14) can be obtained from the following algebraic equations according to the algorithm in Equation (10). Since there is no term about time t or sin ω 1,0 t cos ω 1,0 t after the expansion of Equation (7), Step 2 is ignored.
Step 1: constant term Step 3: the term of cos Step 4: the term of sin 2 ω 1,0 t Step 5: remaining terms → ω 1,1 (t) Given u = 0 is an equilibrium point and the system oscillates within the interval [−A, A], we have Solving Equations (16) and (17) with the conditions Equation (19) yields After the integration of the dynamic frequency ω(t), we obtain According to Equation (13), the frequency of the first order approximate solution is By replacing ω 1,0 with ω 1 in all the nonlinear terms of Equation (21), the first order approximate solution can be expressed as If we stop at the second order dynamic frequency, substituting Equation (18) into Equation (8), we can get the expression of ω 1,2 (t) through a five step balancing process. After the integration of the dynamic frequency ω(t), we obtain By simple calculation, as is illustrated in the process of the first order approximate solution, we can easily get the second order approximate solution which reads The first order approximate period T 1 and the second order approximate period T 2 can be obtained from following relation: According to [5], the exact period of the oscillator is In order to verify the accuracy of the present method, Table 1 shows the comparison of approximate periods T 1 , T 2 , respectively, and the exact period T ex . As illustrated in Table 1, the proposed method lim  (14). Therefore, in the entire range of parameter for a hard spring with γ > 0, the maximal relative error of the period is less than 3.2%. For soft nonlinear oscillators (γ < 0), the approximate period of the present method is still of high accuracy except in a small range of γA 2 near γA 2 = −1, where the heteroclinic orbit with an infinite period appears. For , the exact periodic solution u ex solved by numerical integration of Equation (14), the analytical approximate solutions of different orders u 1 , u 2 calculated by the present method and the approximate solutions of the first two orders u 1 , u 2 obtained by HBM, respectively, as well as their absolute errors, are presented in Figures 1 and 2. These figures show that, the present method provides excellent approximations as compared to the results by the harmonic balance method.

7
In order to verify the accuracy of the present method, Table 1 shows the comparison of approximate periods 1 T , 2 T , respectively, and the exact period ex . T As illustrated in Table 1, the proposed method provides excellent approximations to periods for both small as well as large values of

Oscillator with Coordinate-Dependent Mass
Lev et al. proposed a nonlinear oscillator with coordinate-dependent mass recently [42]. The oscillator is governed by where an overdot represents differentiation with respect to time t and the initial conditions are prescribed as For the above oscillator, it can describe phase transitions in physics, and plays an important role in cosmos-logical model、quark confinement and spinodal. This nonlinear oscillator with a negative coefficient of linear term is difficult to be solved by the perturbation method. In order to apply the dynamic frequency method, Equation (31) can be simplified into the following equivalent form We construct the energy equation of the above equation Submitting Equation (5) into Equation (34), and according to the first order approximation of the present method, we obtain the following equation

Oscillator with Coordinate-Dependent Mass
Lev et al. proposed a nonlinear oscillator with coordinate-dependent mass recently [42]. The oscillator is governed by (1 + αu 2 ) ..
where an overdot represents differentiation with respect to time t and the initial conditions are prescribed as For the above oscillator, it can describe phase transitions in physics, and plays an important role in cosmos-logical model, quark confinement and spinodal. This nonlinear oscillator with a negative coefficient of linear term is difficult to be solved by the perturbation method. In order to apply the dynamic frequency method, Equation (31) can be simplified into the following equivalent form ..
We construct the energy equation of the above equation Submitting Equation (5) into Equation (34), and according to the first order approximation of the present method, we obtain the following equation Here, u = 0 is an equilibrium point and the system oscillates within the interval [−A, A] in the case of α = 0. Thus, we have a 0 = A, b 0 = 0.
Based on the five-step balancing algorithm of the dynamic frequency method, we can get the following algebraic equations according to the algorithm in Equation (10). Since there are no terms about time t and cos ω 1,0 t after the expansion of Equation (35), Step 2 and Step 3 are ignored Step 1: constant term Step 4: the term of sin 2 ω 1,0 t (1 + αA 2 )ω 2 1,0 = (A 2 − 1), Step 5: remaining terms After the integration of ω(t), the first order approximate solution can be obtained Judging from Equation (37), Equations (39) and (40) are valid for the case when It means, when A < 1, there is no periodic solution to Equation (31). In fact, the real number of A could be more or less, but Equation (41) gives an idea of the small amplitude for non-periodic solutions.
According to Yue Wu [43], the homotopy perturbation method (HPM) is an effective technique for such kinds of nonlinear oscillator. In order to compare these two methods, the first approximate period T 1 acquired by using the present method and the first approximate period T calculated by using the homotopy perturbation method are respectively presented in Table 2. The table shows that they have a similar level of accuracy at the range of α near α = 0.1. As the value of α becomes larger, the error between T 1 and T gradually increases. For different values of α, the numerical solution acquired by the Runge-Kutta method and the analytical approximate solutions calculated by the HPM and the present method are respectively presented in Figure 3. The figure shows that the present method provides a superior approximations as compared to numerical solutions.

Modified Energy Balance Method
According to He and other researchers, the solution to nonlinear oscillators can be written as

Modified Energy Balance Method
According to He and other researchers, the solution to nonlinear oscillators can be written as Assuming Equation (42) as the exact solution, the energy equation of the system should be valid for all values of t. However, as Equation (42) is only a type of approximation, the trial Hamilton can only be held at the special collocation point, θ 0 = ω 1,0 t = π/4. That point comes from the balance of kinetic and potential energy at the linear case The so called energy balance method (EBM), is broadly used to simplify the procedure of nonlinear analysis and improve the accuracy of computation. Moreover, some researchers [44,45] have successfully extended this method to non-conservative (damped) oscillatory system, such as Van der Por oscillator. As shown in Equation (43), the initial collocation point θ 0 includes only the linear term. Herein, we use the dynamic frequency to find the more accurate collocation point. Balancing the kinetic and nonlinear potential energy of the Hamilton equation of Equation (1), we can obtain which clearly degenerates to Equation (43) in case M = 1. Substituting the obtained variables (a 0 , b 0 , ω 1,0 ) with ω 1,1 (t) into Equation (44), based on the first order dynamic frequency component ω(t) = ω 1,0 + pω 1,1 (t), it can be numerically solved for the first order collocation point, supposed as θ 1 . That is, obviously superior than the initial point θ 0 by considering the effect of nonlinearities.
Finish the integration in Equation (7), then consider the collocation point θ 1 , and introduce a 1 as a new amplitude variable to replace a 0 . Solving the obtained polynomial for a 1 , and substituting the result into step4 and step5 in Equation (10), we can obtain the modified expressions about ω 1,0 and ω 1,1 (t). Clearly, it avoids the higher order dynamic frequency manipulation in the energy equation so that the computation procedure can be reduced and we can seek a balance between accuracy and workload.
Taking the following Van der Pol equation as an example .. u + ω 2 0 u = (α 2 u 2 + α 3 u 3 + α 4 u 4 + α 5 u 5 ) + (β 0,1 + β 2,1 u 2 ) . u (45) Four different sets of parameters for Equation (45) are listed in Table 3. The new collection point θ 1 is calculated for the purpose of comparison. In Figure 4, we compare the phase diagrams of G 1 to G 4 between the first order dynamic frequency and the modified EBM which transfers the analysis from an isolated linear collocation point to globally effect of nonlinearity. Table 3. Parameter values of the Equation (45).
1,1 equation so that the computation procedure can be reduced and we can seek a balance between accuracy and workload. Taking the following Van der Pol equation as an example Four different sets of parameters for Equation (45) are listed in Table 3. The new collection point 1 θ is calculated for the purpose of comparison. In Figure 4, we compare the phase diagrams of G1 to G4 between the first order dynamic frequency and the modified EBM which transfers the analysis from an isolated linear collocation point to globally effect of nonlinearity.

The Nonlinear Strength Index
In Table 3, a slight difference in the parameter choice will lead to a dramatic change in the nonlinear behavior. Such a change is obviously due to the relationship between the linear and nonlinear parts, which in turn decide the strategies of further computation. Accordingly, one needs an index to classify the strength of nonlinearity so that we can decide to choose the suitable algorithms of computation.
From Equations (43) and (44) In Table 4, we outline a new nonlinear strength index  Table 4.

The Nonlinear Strength Index
In Table 3, a slight difference in the parameter choice will lead to a dramatic change in the nonlinear behavior. Such a change is obviously due to the relationship between the linear and nonlinear parts, which in turn decide the strategies of further computation. Accordingly, one needs an index to classify the strength of nonlinearity so that we can decide to choose the suitable algorithms of computation.
From Equations (43) and (44), the collocation points for the linear and nonlinear Hamilton equations are θ 0 and θ 1 . Suppose θ 1 = π/λ, the absolute value |λ − 4| accordingly represents the effect of nonlinearity, for θ 0 = π/4. In Table 4, we outline a new nonlinear strength index δ = |λ − 4| according to different parameter groups. From G 1 to G 4 , the enlargement of δ synchronizes the intensification of the nonlinearity, as well as the loosing accuracy of the first order analytical approximate solutions. Then, two major nonlinear intervals are logically classified, where the boundary values δ = 0.3 are a coarse estimate, only for the choice of the following analytical strategies. The nonlinear strength indexes and belonging intervals from G 1 to G 4 are presented in Table 4. In Figure 5, we present the flow chart of approximation computation regarding δ. For example, in interval I, it just requires the first order dynamic frequency for the weakly nonlinear feature; in interval II, it demands the second order or other possible kth order dynamic frequency, k ≥ 2, and the modified EBM results improve the accuracy of computation, for its strongly nonlinear character.

Conclusions
In this paper, we present a new dynamic frequency method to solve the periodic solutions for strongly nonlinear oscillators. Based on the analysis of the energy equation, the analytical approximate solutions with high accuracy are obtained through a limited trigonometric balancing progress. Moreover, we set up the strategies to find the higher-order dynamic frequency approximation as well as the modified EBM. Two examples have been solved and the analytical approximate solutions are presented to illustrate the effectiveness and applicability of the present method.
The principal merits of our approach can be summarized as follows: 1. Unlike the harmonic balance method, the influence of the higher harmonic terms is directly embodied in the expression of dynamic frequency which reduces the cumbersome steps of multiple balancing in HB; 2. the new approach is not only valid for strongly nonlinear conservative oscillators but also for nonconservative systems with limit cycles; 3. the lowest order approximations obtained by the present method are actually of high accuracy; 4. the modified energy balance method is introduced to simplify the procredure of high-order analytical computation; and 5. the nonlinear strength index provides possible choices of the avaliable analytical strategies to the general nonlinear systems.
The entire procedure can be regarded as a foundation for further nonlinear analysis, such as for fractional-order nonlinear systems or fluid-structure interaction [46,47]. These will be the topics for further research.

Conclusions
In this paper, we present a new dynamic frequency method to solve the periodic solutions for strongly nonlinear oscillators. Based on the analysis of the energy equation, the analytical approximate solutions with high accuracy are obtained through a limited trigonometric balancing progress. Moreover, we set up the strategies to find the higher-order dynamic frequency approximation as well as the modified EBM. Two examples have been solved and the analytical approximate solutions are presented to illustrate the effectiveness and applicability of the present method.
The principal merits of our approach can be summarized as follows: 1. Unlike the harmonic balance method, the influence of the higher harmonic terms is directly embodied in the expression of dynamic frequency which reduces the cumbersome steps of multiple balancing in HB; 2. the new approach is not only valid for strongly nonlinear conservative oscillators but also for non-conservative systems with limit cycles; 3. the lowest order approximations obtained by the present method are actually of high accuracy; 4. the modified energy balance method is introduced to simplify the procredure of high-order analytical computation; and 5. the nonlinear strength index provides possible choices of the avaliable analytical strategies to the general nonlinear systems.
The entire procedure can be regarded as a foundation for further nonlinear analysis, such as for fractional-order nonlinear systems or fluid-structure interaction [46,47]. These will be the topics for further research.