Homotopy Analysis Method for Nonlinear Dynamical System of an Electrostatically Actuated Microcantilever

Abstract

The homotopy analysis method (HAM) is employed to propose an approach for solving the nonlinear dynamical system of an electrostatically actuated micro-cantilever in MEMS. There are two relative merits of the presented HAM compared with some usual procedures of the HAM. First, a new auxiliary linear operator is constructed. This operator makes it unnecessary to eliminate any secular terms. Furthermore, all the deformation equations are purely linear. Numerical examples show the excellent agreement of the attained solutions with numerical ones. The respective effects of applied voltage, cubic nonlinear stiffness, gap distance, and squeeze film damping on vibration responses are analyzed detailedly.

1. Introduction

Over the last couple of decades, Liao [1–4] described a nonlinear analytical technique which does not require small parameters and thus can be applied to solve nonlinear problems without small or large parameters. This technique is based on homotopy, which is an important part of topology, called the homotopy analysis method (HAM). Its main idea is to construct a class of homotopy in a rather general form by introducing an auxiliary parameter. This parameter can provide us with a convenient way to control the convergence of approximation series and adjust convergence rate and region when necessary. A systematical description of this method was presented in [5]. Liao [5] also studied the convergence properties of the HAM and proved that as long as an HAM series is convergent, it must converge to one solution of the considered problem.

Many microelectro-mechanical systems (MEMSs) inherently contain nonlinearities, such as intrinsic and exterior nonlinearities arising from coupling of different domains [6–8]. Also, there exist mechanical nonlinearities, for example, large deformations, surface contact, creep phenomena, time-dependent masses and nonlinear damping effects, and so forth [9, 10]. It seems fair to say that nonlinear dynamic analysis becomes an increasingly important task in MEMS research and manufacturing.

The nonlinear dynamical behaviors of microcantilever-based instrument in MEMS under various loading conditions have stimulated the curiosities and interests of many researchers [11–15]. For example, the oscillation of an electrostatically actuated microcantilever-based device in MEMS was investigated through a simplified mass-spring-damping model subjected to nonlinear electrostatic force [16–19]. Complex nonlinear terms arising from electrostatic force and from squeeze film damping make it difficult to analyze the system directly using some routine techniques for nonlinear vibrations. For this reason, Zhang and Meng [20] suggested an approximate treatment by expanding the nonlinearities into Taylor series and retaining only the first two terms. The harmonic balance method was then applied to solve the approximate system. Note that the attained results have not been compared with any numerical solutions. The purpose of this study is to seek highly accurate solutions of the aforementioned system on the basis of the HAM.

2. Dynamical System

The electrostatically actuated microcantilever in MEMS as shown in Figure 1 is 4.5𝜇m×80𝜇m×200𝜇m in dimensions [16]. The governing equation describing the motion of the microstructure is [16, 20] 𝑚𝑦^..+𝑐𝑦̇+𝑘𝑦=𝐹_𝐸(𝑡),(2.1) where the superscript denotes the differentiation with respect to time 𝑡,𝑦 the vertical displacement of the microcantilever relative to the origin of the fixed plate, 𝑚 the mass, 𝑘 and 𝑐 the effective spring stiffness and damping coefficient of the simplified system, respectively. According to the parallel plate theory, if the fringe effects at the edges of the plates are ignored [21], the electrostatic force (𝐹_𝐸) generated by applying a voltage 𝑉(𝑡) between the capacitor plates (the fixed plate and the movable plate) can be expressed by𝐹_𝐸=𝜀⁰𝐴2𝑉²(𝑡)(𝑑−𝑦)²,(2.2) where 𝜀₀ is the absolute dielectric constant of vacuum, 𝜀₀=8.5×10⁻¹² N/m, 𝐴 the overlapping area between the two plates, and 𝑑 (μm) is the gap between them. The other parameter values are given as 𝑚=3.5×10⁻¹¹ kg, and 𝑘=0.17 N/m, 𝑐=1.78×10⁻⁶ kg/s and 𝐴=1.6×10⁻⁹ m² [20], while the value of 𝑑 varies.

Figure 1 

A simplified dynamical model of the microcantilever in MEMS.

Introducing the following dimensionless variables𝑥=𝑦𝑑,𝜔⁰=𝑘𝑚,𝑡¹=𝜔₀𝑡,(2.3) then one can rewrite (2.1) and (2.2) as 𝑥^..+𝜁𝑥̇+𝑥=𝑇𝑉²𝑡₁(1−𝑥)²,(2.4) where 𝜁=𝑐/√𝑘𝑚, 𝑇=𝜀₀𝐴/(2𝑘𝑑³), and the superscript (⋅) refers to the differentiation with respect to 𝑡₁.

When the applied voltage, 𝑉(𝑡₁), includes alternating current (ac) voltage, 𝑉₀cos(𝜔𝑡₁), and polarization voltage, 𝑉_𝑃, 𝑉(𝑡₁)=𝑉_𝑃+𝑉₀cos(𝜔𝑡₁) can be considered for analyzing the nonlinearities of the system. In addition, MEMS devices are often characterized by structures that are a few microns in size, separated by micron-sized gaps. At these sizes, air damping dominates over other dissipation. Squeeze film damping may be used to represent the air damping experienced by the moving plates [22]. When considering the effect of squeeze film damping and a cubic nonlinear spring stiffness 𝑘₃𝑥³, system (2.4) is given as 𝑥^..+𝜁+𝜁²(1−𝑥)³𝑥̇+𝑘³𝑥³+𝑥=𝑇𝑉^𝑃+𝑉₀cos𝜔𝑡₁²(1−𝑥)²,(2.5) where 𝜁₂=𝑐₂/(𝑑³√𝑘𝑚) is a nondimensional parameter, and 𝑐₂ is the squeeze film damping (Pa·μm⁴·s).

3. Mathematical Formulation

In this section, the HAM is applied to solve (2.5). First of all, introducing a new time scale 𝜏=𝜔𝑡₁, and multiplying (2.5) by (1−𝑥)³, one can rewrite (2.5) as (1−𝑥)³𝜔²𝑥^+𝜁𝜔𝑥^+𝑥+𝑘₃𝑥³+𝜁₂𝜔𝑥^+𝑇𝑥𝑉_𝑃+𝑉₀cos𝜏²=𝑇𝑉_𝑃+𝑉₀cos𝜏²,(3.1) where the superscript denotes the differentiation with respect to 𝜏. Compared with (2.5), it is very easy to expand (3.1) once 𝑥 is expressed as an HAM series.

3.1. Basic Idea of Homotopy Analysis Method

In the frame of the HAM, the first procedure is to choose an auxiliary linear operator and a nonlinear operator, respectively. A new equation, usually called as the zeroth-order deformation equation, is constructed as (1−𝑝)𝐿[𝑢(𝜏,𝑝)]=ℎ𝑝𝑁[𝑢(𝜏,𝑝)],(3.2) where 𝐿 denotes the chosen auxiliary operator, 𝑁 the nonlinear one, 𝑝∈[0,1] the embedding parameter and ℎ a nonzero auxiliary parameter. The HAM is based on the continuous variation of 𝑢(𝜏,𝑝). When the embedding parameter 𝑝 increases from 0 to 1, 𝑢(𝜏,𝑝) varies from an initial guess to one exact solution of the considered problem 𝑁[𝑢(𝜏,𝑝)]=0.(3.3) According to (3.1), it can be given as𝑁[𝑢(𝜏,𝑝)]=(1−𝑢)³𝜔²𝜕²𝑢𝜕𝜏²+𝜁𝜔𝜕𝑢𝜕𝜏+𝑢+𝑘³𝑢³+𝜁²𝜔𝜕𝑢+𝑇𝑢𝑉_𝑃𝜕𝜏+𝑉₀cos𝜏²−𝑇𝑉_𝑃+𝑉₀cos𝜏².(3.4) When 𝑝=0, (3.2) becomes a linear one, whose exact solution is evident. When 𝑝=1, it is the same as (3.1) provided that 𝑥(𝜏)=𝑢(𝜏,1). Expanding 𝑢(𝜏,𝑝) into a Taylor series 𝑢(𝜏,𝑝)=^∞𝑖=0𝑢_𝑖(𝜏)𝑝^𝑖.(3.5) Note that 𝑢_𝑖(𝜏) is dependent upon ℎ. In fact, the auxiliary parameter ℎ is introduced to make it convenient to control and adjust the convergence of the HAM series (i.e., (3.5)). As long as series (3.5) is convergent at 𝑝=1, one can obtain the 𝑚th-order HAM solution as𝑥(𝜏)=^𝑚𝑖=0𝑢_𝑖(𝜏).(3.6) Substituting (3.4) and (3.5) into (3.2) and equating the coefficients of 𝑝^𝑖 to 0, one obtains𝐿𝑢𝐿𝑢₀=0,(3.7)_𝑛+1−𝑢_𝑛=ℎ𝑅_𝑛𝑢₀,𝑢₁,…,𝑢_𝑛,𝑛=0,1,2,…,(3.8) where𝑅_𝑛=𝜔²^𝑛𝑖=0𝑧_𝑖𝑢^{𝑛−𝑖}+𝜉𝜔^𝑛𝑖=0𝑧_𝑖𝑢^{𝑛−𝑖}+^𝑛𝑖=0𝑧_𝑖𝑢_𝑛−𝑖+𝑘₃^𝑛𝑖=0^{𝑛𝑗=𝑖}^{𝑛𝑙=𝑗}𝑧_𝑖𝑢_𝑗−𝑖𝑢_𝑙−𝑗𝑢_𝑛−𝑙+𝜁₂𝜔𝑢^𝑛+𝑇𝑉_𝑃+𝑉₀cos𝜏²𝑢_𝑛−𝜗_𝑛𝑇𝑉_𝑃+𝑉₀cos𝜏²(3.9) when 𝑛=0, 𝜗_𝑛=1, otherwise 𝜗_𝑛=0; 𝑧_𝑖 corresponds to the expansion of (1−𝑢)³, that is, 𝑧_𝑖=∑^𝑖𝑗=1∑^{𝑖𝑙=𝑗}𝑣_𝑗𝑣_𝑙−𝑗𝑣_𝑖−𝑙, with 𝑣₀=1−𝑢₀ and 𝑣_𝑘=−𝑢_𝑘,𝑘=1,2,…. Equations (3.7) and (3.8) are always linear for every positive integer 𝑛. Importantly, the nonlinear equation has so far been successfully transformed into a series of linear ones. These linear equations can be solved step by step. That is core idea and also the most outstanding advantage of the HAM.

3.2. A Widely Used Linear Operator

One of key procedures of the HAM is to choose an appropriate linear operator. According to Liao [5], it should be chosen on the basis of the so-called rule of solution expression. Periodic/limit cycle solutions can be expressed as Fourier series in 𝜏𝑥=^∞𝑖=0𝛼_𝑖cos(𝑖𝜏)+𝛽_𝑖sin(𝑖𝜏).(3.10) Note that the HAM has been widely used to obtain periodic and/or limit cycle solutions of nonlinear oscillators [23–29]. Particularly, this technique has found applications in microsystems [28] and microbeams subjected applied voltages [29], respectively. To the best of out knowledge, most researchers [23–29] adopted such a following linear operator as𝐿[𝑢(𝜏,𝑝)]=𝜔²𝜕²𝑢(𝜏,𝑝)𝜕𝜏²+𝑢(𝜏,𝑝)(3.11) If (3.11) is employed in this study, (3.7) becomes one homogenous equation in 𝑢₀. This equation has a series of solutions that has to be determined using the first-order deformation equation, that is, (3.8) with 𝑛=0. To this end, an initial guess for 𝑢₀(𝜏) is usually given as𝑢₀=𝑐_0,1cos𝜏+𝑠_0,1sin𝜏,(3.12) where 𝑐_0,1 and 𝑠_0,1 are prior unknown constants.

Due to the rule of solution expression, 𝑢_𝑛 can be given as truncated Fourier series. Substituting 𝑢_𝑛 into 𝑅_𝑛, one has 𝑅_𝑛=^{𝜑(𝑛)𝑖=0}𝐶_𝑛,𝑖cos𝑖𝜏+𝑆_𝑛,𝑖sin𝑖𝜏,(3.13) where 𝜑(𝑛) refers to the highest harmonic. Considering (3.7) and (3.13), cos𝜏 and sin𝜏 can result in the so-called secular terms 𝜏sin𝜏 and 𝜏cos𝜏, respectively. These terms do not comply with the rule of solution expression. One usual procedure to eliminate the secular terms is equating both 𝐶_𝑛,1 and 𝑆_𝑛,1 to 0, that is, 𝐶_𝑛,1=0,𝑆_𝑛,1=0.(3.14) When 𝑛≥1, these equations are all linear. When 𝑛=0, however, (3.14) is nonlinear. Actually, it is essentially the first-order harmonic balance equations [25, 26] governing 𝑐_0,1 and 𝑠_0,1. Well-known, it is a major shortcoming of the harmonic balance method that nonlinear algebraic equations have to be solved.

3.3. A New Linear Operator

Considering that (3.1) contains powers like 𝑥⁶, there will be 𝑐^60,1 and 𝑠^60,1 in (3.14) with 𝑛=0. These nonlinearities can make it rather difficult to determine the initial solution. It is worthy to find another linear auxiliary operator to avoid this obstacle. Consider the following linear operator𝐿[𝑢(𝜏,𝑝)]=𝜔²𝜕²𝑢(𝜏,𝑝)𝜕𝜏²+𝜁+𝜁₂𝜔𝜕𝑢(𝜏,𝑝)𝜕𝜏+𝑢(𝜏,𝑝)−𝑇𝑉^𝑃+𝑉₀cos𝜏²∶=𝐿^__[𝑢(𝜏,𝑝)]−𝑇𝑉_𝑃+𝑉₀cos𝜏².(3.15) Substitution of (3.15) into (3.7) and (3.8) results into𝐿^__𝑢₀(𝜏)−𝑇𝑉_𝑃+𝑉₀cos𝜏²=0,(3.16)_{__}𝐿𝑢_𝑛+1−𝜎_𝑛𝑢_𝑛=ℎ𝑅_𝑛𝑢₀,𝑢₁,…,𝑢_𝑛,𝑛=0,1,2,…,(3.17) where 𝜎₀=0, otherwise 𝜎_𝑛=1 when 𝑛≥1. Equations (3.16) and (3.17) provide us with a purely linear algorithm, without an additional task of eliminating the secular terms as long as 𝜁+𝜁₂≠0. Since (3.16) is an inhomogeneous equation whose exact solution can be determined analytically, there is also no need to choose an initial guess for 𝑢₀.

3.4. A Means for Saving Computational Effort

The number of harmonics in the HAM approximations increases acceleratedly, because of the existence of the powers like 𝑥⁶ in (3.1). Denote the number of harmonics in 𝑢_𝑛 as Γ_𝑛. The starting solution (𝑢₀) provided by (3.16) contains the first two harmonics, that is, Γ₀=2. In the equation governing 𝑢₁, the highest powers is 𝑢⁶⁰, hence one knows that Γ₁=12. As 𝑛 increases further, Γ₂=22 and Γ₃=32 can be obtained. In most cases, periodic or limit cycle oscillations are dominated by several lower harmonics. High harmonics contribute to the improvement of the precision quite limitedly, compared with the lower ones. Therefore, it seems to be unworthy to improve the accuracy slightly with such a vast increasing of computational cost.

Therefore, the high-order HAM approximations are obtained with a given number (𝑁) of harmonics, while neglecting all harmonics higher than the 𝑁th one. Denote the 𝑛th-order HAM approximation as 𝑢_𝑛=^𝑀𝑖=0𝑐_𝑛.𝑖cos(𝑖𝜏)+𝑠_𝑛,𝑖sin(𝑖𝜏),𝑛=0,1,2,…,(3.18) where 𝑀 increases with 𝑛 increasing. Once 𝑀>𝑁, retain the first to the 𝑁th harmonics in 𝑢_𝑛, and rewrite it as𝑢_𝑛=^𝑁𝑖=0𝑐_𝑛,𝑖cos(𝑖𝜏)+𝑠_𝑛.𝑖sin(𝑖𝜏),𝑛=0,1,2,….(3.19) When seeking higher-order HAM approximations, (3.19) but not (3.18) is employed. Through this means, an HAM solution with 𝑁 harmonics can be obtained. Interestingly, as long as this solution converges, it must converge to one harmonic balance solution, as revealed in [30].

4. Results and Discussion

Theoretically, analytical solutions of 𝑢_𝑖 can be obtained by solving (3.16) and (3.17) step by step, since all equations are linear. However, their expressions can rapidly become complicated as 𝑛 increasing. It is rather easy to obtain semi-analytical solutions when all parameters are given. The HAM solutions are compared with numerical results. All the numerical solutions are obtained via integrating (2.5) using the RK method, with the initial values as 𝑥(0)=0 and 𝑥̇(0)=0.

Let 𝑇_𝑁,𝑛 denote the CPU time needed to obtain the 𝑛th-order HAM approximation (𝑢_𝑛) with 𝑁 harmonics, and 𝑇_𝑛 is the counterpart when no harmonic is neglected. In addition, let 𝐴_𝑥 be the maximum value (amplitude) of the 𝑛th-order HAM solution (i.e., ∑^𝑛𝑖=0𝑢_𝑖) with 𝑁 harmonics retained, and 𝐵_𝑥 the counterpart. From Figure 2, one can see the convergence of 𝐴_𝑥 and 𝐵_𝑥 to the numerical one as 𝑛 increasing. It implies that the accuracy will not be lost by retaining only the first 10 harmonics (𝑁=10). Also shown in Figure 2 is the ratio between 𝑇_𝑁,𝑛 and 𝑇_𝑛. It rises rapidly when 𝑛 increases. Therefore, a lot of computational effort can be saved without losing substantial accuracy, through neglecting some high harmonics. Note that, all the following HAM solutions are obtained by employing (3.19) rather than (3.18).

(a)

(b)

Figure 2 

(a) the HAM solutions; (b) the ratio between 𝑇𝑁,𝑛 and 𝑇𝑛, 𝑁=10, and ℎ=−1. Parameter values are given as 𝑉𝑃=1, 𝑉0=2, 𝑑0=1, 𝜔=0.5, 𝑘3=1, and 𝑐2=1.

A preliminary procedure to obtain an HAM solution consists in the choosing of the auxiliary parameter ℎ. Liao [23] suggested an ℎ-curve approach for choosing a proper value of ℎ, which can guarantee both the convergence of HAM series and a relatively rapid convergence rate. Figure 3 shows the relative errors versus varying ℎ. These errors refer to the percentage-wise differences between the amplitudes (𝐴_𝑥) of the 20th-order HAM solutions with 10 harmonics and the numerical results. According to the figure, an appropriate range for ℎ can be chosen from −2 to −0.5 or so. When the absolute value of ℎ is too large, for instance more than 2, it can probably make the HAM series divergent. On the other hand, as long as ℎ approaches 0 enough, the HAM series converge by a relatively slow convergence rate.

Figure 3 

Relative errors versus ℎ between 20th-order HAM solutions and numerical ones. Parameter values are 𝑐2=2, 𝜔=1, 𝑉𝑝=1, and 𝑑=2, respectively.

Figure 4 shows a phase plane of system (2.4). The solutions provided by the HAM approach the numerical one when 𝑚 increases. The 20th-order HAM solution is nearly the same as the numerical one. In fact, solutions can be obtained to any desired accuracy without additional difficulty.

Figure 4 

Comparison between numerical solutions and HAM results ℎ=−1, where heavy dots denote the former and solid lines the latter. Parameter values are given as 𝑉𝑃=1, 𝑉0=3, 𝑑=1.2, 𝑘3=𝑐2=0, and 𝜔=1.

When a nonlinear spring and/or squeezefilm damping are taken into account, the proposed algorithm can still provide very accurate solutions, as shown in Figures 5 and 6. Both the figures show clearly the convergence of the HAM solutions to the numerical ones. It is worthy of emphasizing that, the vibrations described by these phase planes are nonharmonic, namely the solutions cannot be approximated by the first harmonic alone. That is probably because high harmonics contribute to the responses to an nonnegligible extent. The phase plane shown in Figure 6 even has a loop. Even though, the presented HAM can still track them with excellent precision. These illustrative examples demonstrate the feasibility and efficiency of the presented approach in analyzing nonlinear dynamical behaviors.

Figure 5 

Comparison between numerical solutions and HAM solutions ℎ=−1, where heavy dots denote the former and solid lines the latter. Parameter values are given as 𝑉𝑃=1, 𝑉0=5, 𝑑=1.5, 𝑘3=0.5, 𝑐2=0, and 𝜔=1.

Figure 6 

Comparison between numerical solutions and HAM solutions ℎ=−1, where heavy dots denote the former and solid lines the latter. Parameter values are given as 𝑉𝑃=0.5, 𝑉0=2, 𝑑=1, 𝑘3=0.5, 𝑐2=2, and 𝜔=0.5.

Frequency-response curve is of substantial importance for dynamic analysis, when external excitations are included. Figures 7, 8, and 9 show the curves of the amplitudes of periodic solutions versus varying 𝜔, for different values of 𝑉_𝑃, 𝑉₀ and 𝑑, respectively. The 20th-order HAM solutions are in excellent agreement with the numerical ones. The amplitudes decrease with 𝜔 increasing for all the three cases. According to (2.5), 𝑉_𝑃 and 𝑉₀ correspond to the magnitude of the external excitation, respectively. Also, 𝑑 is connected with the external excitation through 𝑇=𝜀₀𝐴/(2𝑘𝑑³), which implies that 𝑇 increases when 𝑑 decreases. That is why the amplitude increases when 𝑉_𝑃(𝑉₀) increases or 𝑑 decreases, as shown in these figures.

Figure 7 

The vibrations amplitudes versus 𝜔, where heavy dots denote numerical results and solid lines the 20th-order HAM solutions with 𝑁=10, ℎ=−2, 𝑉0=2, 𝑑0=1, 𝑘3=0.5, and 𝑐2=2.

Figure 8 

The solutions amplitudes versus 𝜔, where heavy dots denote numerical results and solid lines the 20th-order HAM solutions with 𝑁=10, ℎ=−2, 𝑉𝑃=0.5, 𝑑0=1, 𝑘3=0.5, and 𝑐2=2.

Figure 9 

The solutions amplitudes versus 𝜔, where heavy dots denote numerical results and solid lines the 20th-order HAM solutions with 𝑁=10, ℎ=−2, 𝑉𝑃=1, 𝑉0=2, 𝑐2=2, and 𝑘3=0.5.

Figures 10 and 11 are presented to investigate the respective effects on the periodic responses of the cubic stiffness coefficient and the squeeze film damping. Nice agreement of the 20th-order HAM solution to the numerical results can also be observed. As they show, respectively, the larger 𝑘₃ or 𝑐₂ is, the smaller the amplitudes become. According to (2.5), 𝑘₃𝑥³ is the nonlinear restoring force, hence 𝑘₃𝑥² can be considered as an equivalent linear stiffness coefficient. Its value keeps nonnegative though varying. The increasing of 𝑘₃ implies an increasing of the stiffness, so that the magnitude of response decreases. As concerning 𝑐₂, it denotes the coefficient of a positive nonlinear damping. Generally speaking, positive damping can suppress the vibration to some extent.

Figure 10 

The solutions amplitudes versus 𝜔, where heavy dots denote numerical results and solid lines the 20th-order HAM solutions with 𝑁=10, ℎ=−2, 𝑉𝑃=1, 𝑉0=2, 𝑑0=1, and 𝑐2=2.

Figure 11 

The solutions amplitudes versus 𝜔 where heavy dots denote numerical results and solid lines the 20th-order HAM solutions with 𝑁=10, ℎ=−2, 𝑉𝑃=1, 𝑉0=2, 𝑑0=1, and 𝑘3=0.5.

It is worthy to point out, the influence of squeeze film damping on vibration amplitudes reduces rapidly when 𝜔 increases. When 𝜔 is relative large (𝜔>1.5), the decreasing of 𝑐₂ makes a very little difference of the amplitudes.

5. Conclusions and Remarks

Based on the homotopy analysis method, we have proposed an approach to solve the nonlinear dynamical system of an electrostatically actuated microcantilever. This approach was further improved by restricting the number of harmonics as a given number 𝑁. It turns out that the improved approach can save a large amount of computational effort without reducing accuracy. Numerical examples validate the proposed approach. Approximations for periodic solutions are obtained very precisely. Using the presented approach, the frequency-response curves are obtained and discussed in details. The feasibility and efficiency, illustrated by the numerical examples, imply that we could expect the presented approach be applicable in more nonlinear vibration problems, especially those arising in MEMS.

As mentioned above, we have presented an auxiliary operator (3.15) that is different from another usually-adopted one (3.11). It is advantageous to use (3.15) rather than (3.11) for the considered problem. First, there is no need to choose an initial guess for the zeroth-order HAM approximation if one employs linear operator (3.15). Second, if (3.11) is adopted, a series of algebraic equations have to be solved from time to time to eliminate secular terms when seeking every order HAM approximation. In many occasions, the first equation is nonlinear. In our procedures, there is no requirement of doing so. Last but not least, the proposed procedures are purely linear.

Acknowledgments

This paper is supported by the National Natural Science Foundation of China (10102023, 10972241, 11002088 and 10732060), Doctoral Program Foundation of Ministry of Education of China (20090171110042), China Postdoctoral Science Foundation (20090460628), and Shanghai Postdoctoral Science Foundation (10R21413600). Project supported by the Research Fund of State Key Lab of MSV, China (Grant no. MSV-MS-2010-9).