An Iteration Scheme Suitable for Solving Limit Cycles of Nonsmooth Dynamical Systems

Abstract

The Mickens iteration method (MIM) is modified to solve self-excited systems containing nonsmooth nonlinearities and/or nonlinear damping terms. If the MIM is implemented routinely, the unknown frequency and amplitude of limit cycle (LC) would couple to each other in complicated nonlinear algebraic equations at each iteration. It is cumbersome to solve these algebraic equations, especially for nonsmooth systems. In the modified procedures, the unknown frequency is substituted by the determined value obtained at the previous iteration. By this means, the frequency is decoupled from the nonlinear terms. Numerical examples show that the LCs obtained by the modified MIM agree well with numerical results. The presented method is very suitable for solving self-excited systems, especially those with nonlinear damping and nonsmooth nonlinearities.

1. Introduction

Recent years have witnessed the wide applications of iteration techniques, such as the Mickens iteration method (MIM) [1–6] and the variational iteration method [7–10]. In order to improve the efficiency, Lim and Wu [11], Marinca and Herisanu [12], and Hu [13] modified the MIM, respectively.

In principle, the approximations can be obtained to any desired accuracy by the MIM as long as the iteration proceeds. In the MIM, algebraic equations are introduced at each iteration to eliminate the so-called secular terms. In the applications in conservative oscillators, the algebraic equations are linear. Nonlinear damping and nonsmooth terms appear widely in many dynamical systems [14, 15]. As for the oscillators with nonlinear damping terms, however, very complicated nonlinear algebraic equations have to be solved at each iteration [16]. Moreover, the algebraic equations cannot be deduced for systems with nonsmooth nonlinearities. It is necessary and worthwhile, therefore, to propose some approaches to simplify the MIM. This paper will present a modified iteration algorithm by decoupling the unknown frequency from nonlinear terms.

2. A Modified MIM

Consider a self-excited oscillator 𝑥^..+𝑓(𝑥,𝑥̇)=0, (1) where the superscript denotes the differentiation with respect to time 𝑡 and 𝑓(𝑥,𝑥̇) is a nonlinear term with damping terms. Assume that system (1) has at least one limit cycle (LC) solution. Since the LC frequency and amplitude are independent of initial conditions, they should be considered as unknowns to be determined at every iteration. Denote the angular frequency as 𝜔 and introduce the transformation as 𝜏=𝜔𝑡; thus, we rewrite (1) as 𝜔²𝑥^′′+𝜔²𝑥=𝜔²𝑥−𝑓𝑥,𝜔𝑥^′ (2) subject to the following initial conditions: 𝑥(0)=𝛼,𝑥^′(0)=0, (3) where 𝛼 is the unknown LC amplitude and the superscript denotes the differentiation with respect to 𝜏. Note that 𝛼 will be approximated as a series {𝛼_𝑘} by eliminating the secular terms at each iteration stage. In order to obtain the LC, the MIM [16] can be given as 𝜔²_𝑘−1𝑥^′′_𝑘+𝑥_𝑘=𝜔²_𝑘−1𝑥_𝑘−1−𝑓𝑥_𝑘−1,𝜔_𝑘−1𝑥^′_𝑘−1,𝑘=1,2,…, (4) with the initial conditions being rewritten at each iteration as 𝑥_𝑘(0)=𝛼_𝑘,𝑥^′_𝑘(0)=0. (5) Note that the coefficient of the first harmonic in 𝑥_𝑘−1 remains still to be an unknown, that is, 𝛼_𝑘−1. This unknown will couple with the unknown frequency, 𝜔_𝑘−1, which will result in a coupled nonlinear term (i.e., 𝜔_𝑘−1𝑥^′_𝑘−1) in the right side of (4). If higher powers of 𝑥^′ exist, these terms will lead to very complicated functions in 𝜔_𝑘−1 and 𝛼_𝑘−1. Different from conservative systems, 𝜔²_𝑘−1 can no longer be considered as an independent unknown. In order to simplify the MIM, therefore, a modified scheme is proposed as 𝜔²_𝑘−1𝑥^′′_𝑘+𝑥_𝑘=𝜔²_𝑘−1𝑥_𝑘−1−𝑓𝑥_𝑘−1,𝜔_𝑘−2𝑥^′_𝑘−1,𝑘=1,2,…. (6) As 𝑘=1, we choose 𝜔₋₁=𝜔₀. In the 𝑘th iteration, 𝜔_𝑘−2 is a given constant that is obtained at the (𝑘−1)th iteration. The square of the unknown frequency, that is, 𝜔²_𝑘−1, can be treated as an independent parameter, because 𝜔_𝑘−1 appears only in 𝜔²_𝑘−1𝑥_𝑘−1. According to the initial conditions, the starting iteration solution can be chosen as 𝑥₀(𝜏)=𝛼₀cos𝜏. (7) It is obvious that as long as the series {𝜔_𝑘,𝑘=1,2,…} and {𝑥_𝑘,𝑘=1,2,…} are convergent, they must converge to the exact solutions. The right-hand side of (6) can be expressed by Fourier series as 𝜔²_𝑘−1𝑥_𝑘−1−𝑓𝑥_𝑘−1,𝜔_𝑘−2𝑥^′_𝑘−1=^𝜑(𝑘)_𝑖=1𝑐_{𝑘−1,𝑖}𝜔²_𝑘−1,𝛼_𝑘−1cos(𝑖𝜏)+𝑠_{𝑘−1,𝑖}𝜔²_𝑘−1,𝛼_𝑘−1sin(𝑖𝜏), (8) where the harmonic coefficients 𝑐_{𝑘−1,𝑖}(𝜔²_𝑘−1,𝛼_𝑘−1) and 𝑠_{𝑘−1,𝑖}(𝜔²_𝑘−1,𝛼_𝑘−1) are functions in 𝜔²_𝑘−1 and 𝛼_𝑘−1. Here, 𝜑(𝑘) is a positive integer denoting the order of the highest harmonic. Approximations 𝜔_𝑘−1 and 𝛼_𝑘−1 are determined by eliminating the secular terms, that is, letting 𝑐_𝑘−1,1𝜔²_𝑘−1,𝛼_𝑘−1=0,𝑠_𝑘−1,1𝜔²_𝑘−1,𝛼_𝑘−1=0,𝑘=1,2,…. (9) These equations can be solved analytically if 𝜔²_𝑘−1𝛼_𝑘−1 is considered as an independent unknown. They can also be numerically solved by Newton-Raphson method. The latter is employed in this study.

Different from the existing procedures [16], as 𝑘 increases, 𝜔²_𝑘−1 is always an independent unknown in the modified MIM. Moreover, unknown 𝜔_𝑘−1 does not couple with nonlinear terms. It simplifies the MIM to a large extent, as shown later.

3. Numerical Examples

Example 1 (system with nonlinear damping terms). The van der Pol equation is chosen to illustrate the previous procedures more clearly: 𝑥^..+𝑥+𝜀𝑥²−1𝑥̇=0, (10) where 𝜀 is a given constant. As known, (10) has a stable LC solution when 𝜀>0 while an unstable one when 𝜀<0.
According to the modified MIM, the corresponding iteration scheme is given as 𝑥^.._𝑘+𝜔²_𝑘−1𝑥_𝑘=𝜔²_𝑘−1𝑥_𝑘−1−𝑥_𝑘−1−𝜀𝑥²_𝑘−1−1𝑥̇_𝑘,𝑘=1,2,…. (11) Introducing a new time variable 𝜏=𝜔_𝑘−1𝑡 at each iteration stage, we rewrite (11) as 𝜔²_𝑘−1𝑥^′′_𝑘+𝑥_𝑘=𝜔²_𝑘−1𝑥_𝑘−1−𝑥_𝑘−1−𝜀𝜔_𝑘−1𝑥²_𝑘−1−1𝑥^′_𝑘−1, (12) where the superscript denotes the derivative with respect to 𝜏. The iteration algorithm begins with an initial solution 𝑥₀(𝜏)=𝛼₀cos𝜏. (13) Then, we obtain the governing equations in 𝑥₁(𝜏) as 𝜔²₀𝑥^′′₁+𝑥₁=𝜔²₀𝛼₀−𝛼₀cos𝜏−𝜀𝛼₀𝜔₀𝛼−𝜀³⁰𝜔₀4𝛼sin𝜏+𝜀³⁰𝜔₀4sin3𝜏,𝑥₁(0)=𝛼₁,𝑥^′₁(0)=0. (14) Equating the coefficients of cos𝜏 and sin𝜏 to zeros results into 𝜔²₀𝛼₀−𝛼₀=0,𝜀𝛼₀𝜔₀𝛼−𝜀³⁰𝜔₀4=0 (15) which yields that 𝜔₀=1 and 𝛼₀=2. Substituting them into (16), we have 𝑥^′′₁+𝑥₁=2𝜀sin3𝜏. (16) Considering initial conditions (14), we can obtain 𝑥₁=𝛼₁cos𝜏+3𝜀sin𝜏4−𝜀sin3𝜏4, (17) where 𝛼₁ is to be determined at the next iteration stage.
According to iterative scheme (11), the equation in 𝑥₂(𝜏) is deduced as 𝜔²₁𝑥^′′₂+𝑥₂=⁹_𝑖=1𝑐_1,𝑖𝜔²₁,𝛼₁cos𝑖𝜔₁𝑡+𝑠_1,𝑖𝜔²₁,𝛼₁sin𝑖𝜔₁𝑡. (18) Equate the coefficients of cos𝜏 and sin𝜏 to zeros: 3𝜔²₁4−15𝛼₁−𝛼644³¹−34𝛼=0,₁𝜔²₁−𝛼₁−𝛼³₁8+75128=0. (19) By solving (19) numerically, we can determine 𝛼₁ and 𝜔²₁. Here, we obtain the second-order approximation and expand it as 𝜔₂^IS𝜀=1−²+1629𝜀⁴𝜀2048+𝑜⁴. (20) According to [17], the Lindstedt-Poincare (LP) method provides the second- and forth-order approximate frequency 𝜔₂^LP=1−𝜀²/16 and 𝜔₄^LP=1−𝜀²/16+17𝜀⁴/3072, respectively. The attained approximation agrees well with the 4th-order LP solution. Table 1 indicates that 𝜔₂^IS is more accurate than 𝜔₂^LP when compared with 𝜔₄^LP.

Table 1 

Comparison of the second-order frequency obtained by IS and LP method with the forth-order approximation obtained by LP method, when 𝜀=1.

Figure 1 shows the comparison of the phase planes between iteration solutions (𝑥_𝑘) and numerical result. Rapid convergence of 𝑥_𝑘 to the numerical result can be observed. Note that all numerical solutions are obtained by the fourth-order Runge-Kutta (RK) integration method. When |𝜀|>1, the iteration procedure presented by Chen and Liu [16] does not converge. This is probably the difference between the starting function (𝑥₀=𝛼₀cos𝑡) and the exact solution is too large. The modified MIM is still effective for |𝜀|≤1.5. As Figure 2 shows, the LC solution with 𝜀=1.5 obtained by the presented method is in excellent agreement with numerical one. It is necessary to point out that the presented method is able to track unstable LCs, whereas the RK method is not.

Figure 1 

The LC solutions of system (10) with 𝜀=1 obtained by the modified MIM and RK method, respectively. The 𝑘th-order solution (𝑥𝑘) is represented by IS𝑘.

Figure 2 

The LC solutions of system (10) with 𝜀=1.5 obtained by the modified MIM, RK method, and LP method, respectively. The iteration solution is denoted as IS and the 𝑘th-order LP approximation as LP𝑘.

Also plotted in Figure 2 are the results provided by the LP method [17]. The iteration results are much more precise than the 2nd and 4th-order LP approximations.

It is necessary to point out that the presented method is able to track unstable LCs, whereas the RK method is not. Figure 3 shows an unstable LC of the van der Pol equation with 𝜀=−1 obtained by the presented method. As shown, the RK begins at the LC; however, the solution curve converges to the equilibrium.

Figure 3 

Comparison of the LC solution of system (10) with 𝜀=−1, provided by the modified MIM, and by RK method, respectively.

In order to further demonstrate the merit of the modified MIM when applied to problems with nonlinear damping terms, we consider the following self-excited system [18]: 𝑥^..+𝑥+𝜀(𝑥̇)²−1𝑥̇+𝑓(𝑥̇)=0. (21) The nonlinear term contains high powers of 𝑥̇, that is, 𝑓(𝑥̇)=(𝑥̇)³. If the original MIM is employed, the algebraic equations governing 𝜔 will become very complicated. Therefore, it is necessary to employ the modified approach. Figure 4 indicates that the approximations obtained by the presented method converge rapidly to the numerical solution as 𝑘 increased.

Figure 4 

The LC solutions of system (21) with 𝜀=1 obtained by the modified MIM and RK method, respectively. The 𝑘th-order solution (𝑥𝑘) is represented by IS𝑘.

Example 2 (system with nonsmooth nonlinearity). The modified MVIM is further applied to nonsmooth dynamical system expressed as 𝑥^..+𝑓(𝑥,𝑥̇)+𝜂𝑔(𝑥)=0. (22) Here, 𝑓(𝑥,𝑥̇) is a nonlinear damping term, and 𝑔(𝑥) is a nonsmooth function. If substituting the 𝑥_𝑘−1 into 𝑔(𝑥), on account of 𝑥_𝑘−1 contained unknown quantities (𝛼_𝑘−1), so 𝑔(𝑥) can not be expanded as Fourier progression by numerical integration. To this end, (9) cannot be deduced by eliminating the secular terms. Likewise, we present the following iteration scheme: 𝜔²_𝑘−1𝑥̈_𝑘+𝑥_𝑘=𝜔²_𝑘−1𝑥_𝑘−1+𝑓𝑥_𝑘−1,𝜔_𝑘−2𝑥̇_𝑘−1+𝜂𝑔𝑥_𝑘−2. (23) In this scheme, 𝑔(𝑥_𝑘−2) can be expanded as a Fourier series since 𝛼_𝑘−2 has been determined at the previous iteration.
Let us consider a van der Pol type oscillator with a nonsmooth function as 𝑥̈+𝑥+𝜀1−𝑥²𝑥̇+𝜂𝑔(𝑥)=0 (24) with 𝑔(𝑥)=𝑥−1𝑥≥10−1<𝑥<1𝑥+1𝑥≤−1. (25) Figure 5 shows the LC of system (14) with 𝜀=1 and 𝜂=0.5. The 5th-order approximations obtained by the presented method agree well with the numerical solution when the nonsmooth term is given as 𝑔(𝑥)=sgn(𝑥)=1,𝑥>00,𝑥=0−1,𝑥<0. (26) The LC can also be obtained very accurate, as Figure 6 shows.

Figure 5 

LC solutions of system (24) with (25) (𝜀=1, 𝜂=0.5) provided by the modified MIM and by RK method, respectively.

Figure 6 

LC solutions of system (24) with (26) (𝜀=1, 𝜂=0.5) provided by the modified MIM and by RK method, respectively.

4. Conclusions

The Mickens iteration method (MIM) has been modified, so that it is suitable for solving LC solutions of self-excited systems with nonsmooth and/or damping nonlinearities. Different from the routinely-used MIM, the modified method decouples the unknown frequency from nonlinear terms. This modification simplifies the MIM significantly. Numerical examples show the feasibility and validity of the presented method, which implies that it could be applicable to more nonlinear dynamical systems, especially those with nonlinear damping terms and nonsmooth nonlinearities.

Acknowledgment

This work is supported by the National Natural Science Foundation of China (11002088, 11272361, 11172333), Doctoral Program Foundation of Ministry of Education of China (20130171110039), Guangdong Province Natural Science Foundation (S2012040007920, S2013010013802), Fundamental Research Funds for the Central Universities (13lgzd06), and the Guangdong Province Science and Technology Program (2012A030200011).