A coupling successive approximation method for solving duffing equation and its application

Volume 36 Number 2 2 2014 Vietnam Journal of Mechanics, VAST, Vol. 36, No. 2 (2014), pp. 77 – 93 A COUPLING SUCCESSIVE APPROXIMATION METHOD FOR SOLVING DUFFING EQUATION AND ITS APPLICATION Dao Huy Bich1, Nguyen Dang Bich2,∗ 1Hanoi University of Science, VNU, Vietnam 2Institute for Building Science and Technology (IBST), Hanoi, Vietnam ∗E-mail: dangbichnguyen@gmail.com Received March 15, 2014 Abstract. The paper proposes an algorithm to solve a general Duffing equation, in which a process of transforming the initial equation to a resulting equation is proposed, and then the coupling successive approximation method is applied to solve the resulting equation. By using this algorithm a special physical factor and complex-valued solutions to the general Duffing equation are revealed. The proposed algorithm does not use any assumption of small parameters in the equation solving. The coupling successive proce- dure provides an analytic approximated solution in both real-valued or complex-valued solution. The procedure also reveals a formula to evaluate the vibration frequency, ϕ, of the non-linear equation. Since the first approximation solution is in a closed-form, the chaos index of the general Duffing equation and the chaotic characteristics of solutions can be predicted. Some examples are used to illustrate the proposed method. In the case of chaotic solution, the Pointcaré conjecture is used for solution verification. Keywords : General Duffing equation, coupling successive approximation method, chaos index, chaotic structures of solutions. 1. INTRODUCTION Many dynamical problems lead to a general Duffing equation [1–6] with not only the third order but also second order nonlinearities. Recently, there have been increasing researches dealing with the Duffing equation with a set of chaotic solutions. The commonly index used to recognize the chaotic solutions of the Duffing equation is Liapunov index [7], which requires a complex calculation. Therefore, there is a need to develop another index using an easier calculation. There are many available methods to solve nonlinear differential equations in gen- eral, and the Duffing equation in particular. For examples: - The Rungen Kutta [8] and Newmark [9] methods are two numerical methods, which require no additional assumptions in equation solved. Applying these methods, some authors were able to find and assess chaotic solutions to the Duffing equation. However, 78 Dao Huy Bich, Nguyen Dang Bich these solutions may not converge [10] and there are difficulties in assessing properties of parameter-dependent solutions. - The perturbation method [11], the asymptotic method [12], the averaging method [13] and the extended Galerkin [14] method are four approximate analytical methods that requires an assumption of small parameters in applying the algorithm to the equation and convergence. However, in solving mechanical and dynamical problems [1–6], coefficients of nonlinear terms and coefficients of exciting forces are rather large. The assumption of small parameters also leads to another limitation in the four methods, where only periodic solutions could be found and examined, but not chaotic solutions. The methods used in Refs. [8–14] were able to find and assess the real-valued solu- tions of the Duffing equation. These findings are considered sufficient for linear differential equations since the linear combination of the real and imaginary components of the so- lutions is also a solution to this type of equations. However, for nonlinear differential equations, these findings are not sufficient as the linear combination is not a solution to these equations [15]. Moreover, for the complex-valued solutions of nonlinear differential equations, there exists a phase space instead of phase plane as for the real-valued solutions [15]. In order to overcome the shortcoming, many new techniques appeared in the lit- erature such as: variational iteration method [16, 17], energy balance method [18, 19], parameter-expansion method [20], Hamiltonian approach [21, 22], variational approach [23, 24], homotopy analysis method [25], . . . When most of these methods are applied to free nonlinear vibration equation the convergence of method is still open. They do not indicate the way to find complex-valued solutions of the equation if that exist. This paper proposes an algorithm and an iterative process to find the real-valued solution or complex-valued solution of the general Duffing equation. The focus of interests in the paper includes: Finding an algorithm to transform the initial equation to a resulting equation; proposing a coupling successive approximation method to solve the resulting equation; finding an analytical approximation solution in real-valued or complex-valued solutions; building the formulae for calculating the oscillation frequency of a system of non-linear equations, and providing a simple indication of chaotic solutions to the Duffing equation; describing chaotic solutions in an analytical form, examining the characteristics of chaotic solutions and using the Poincaré conjecture for verifying the chaotic properties of the so- lutions. 2. THE ALGORITHM TO TRANSFORM THE INITIAL EQUATION TO THE RESULTING EQUATION 2.1. The algorithm to transform the initial equation Consider a general Duffing equation as x¨+ 2νx˙+ λx3 + 2qx2 + kx = p cosω t. (1) A coupling successive approximation method for solving Duffing equation . . . 79 Adding a fixed parameter σ to both sides of Eq. (1), then dividing both sides by x + d, we have 1 x+ d ( x¨+ 2νx˙+ λx3 + 2qx2 + kx+ σ ) = σ + p cosω t x+ d , (2) where d is a fixed parameter, x+ d is assumed to be different from zero and checked after x is found. Eq. (2) can be rewritten as d dt ( u˙ u ) + u˙2 u2 + d dt ( 1 2 υ˙ υ ) − ( 1 2 υ˙ υ )2 = σ + p cosω t x+ d , (3) in which u˙ u = a1 x˙ x+ d + a2x+ a3 + b x+ d , (4) 1 2 υ˙ υ = b1 x˙ x+ d + b2x+ b3 + b x+ d (5) and with the conditions that parameters σ, d, a1, a2, a3, b1, b2, b3, b are suitably determined. Applying Eqs. (4) and (5) into Eq. (3), then combining the equivalent equation with Eq. (2), we obtain a system of algebraic nonlinear equations of parameters a1 − b1 = 1, a2 + b2 + 2a1a2 − 2b1b2 = 0, (a2 + b2) d+ 2a1a3 − 2b1b3 = 2ν (a1 + b1) , a22 − b 2 2 = λ (a1 + b1) ,( a22 − b 2 2 ) d+ 2a2a3 − 2b2b3 = 2q (a1 + b1) , (2a2a3 − 2b2b3) d+ a 3 3 − b 3 3 + (2a2 − 2b2) b = k (a1 + b1) ,( a33 − b 3 3 ) d+ (2a3 − 2b3) b = σ (a1 + b1) . (6) Solving the system of Eqs. (6), the parameters are given by b1 = a1 − 1, a2 = (2a1 − 3) √ − λ 8 , b2 = (2a1 + 1) √ − λ 8 , a3 = 1 3 √ −λ8 [ (2a1 + 1) ν √ − λ 8 − (a1 − 1) q + 3 8 (2a1 − 1) λd ] , b3 = 1 3 √ −λ8 [ (2a1 − 3) ν √ − λ 8 − a1q + 3 8 (2a1 − 1)λd ] , b = 2a1 − 1 8 √ −λ8 [ (2q − λd) d+ 8 9λ ( q + 4ν √ − λ 8 )( q − 2ν √ − λ 8 ) − 2 3 ( q+4ν √ − λ 8 ) d− k ] , σ = − 2 3λ ( q + 4ν √ − λ 8 )[ 8 9λ ( q + 4ν √ − λ 8 )( q − 2ν √ − λ 8 ) − k ] . (7) 80 Dao Huy Bich, Nguyen Dang Bich Given Eq. (4) and Eq. (5) and the parameters specified by Eq. (7) one can see that Eq. (3) and Eq. (1) are totally equivalent. Note that, there are still two arbitrary parameters, a1 and d in Eq. (7). For simplicity, a1 and d are suitably selected. First, d is selected so that b = 0, i.e. d is selected from the equation (2q − λd) d+ 8 9λ ( q + 4ν √ − λ 8 )( q − 2ν √ − λ 8 ) − 2 3 ( q + 4ν √ − λ 8 ) d− k = 0. Solving the equation yields d = 2 3λ ( q − 2ν √ − λ 8 ) ± θ 2 √ −λ8 , (8) where θ = [ 1 2 ( k − 4 3 q2 λ − 2 3 ν2 )] 1 2 . (9) Next, a1 is selected so that b1 = 0. Thus, a1 = 1, a2 = − √ − λ 8 , b2 = 3 √ − λ 8 , a3 = 5 6 ν ∓ θ 2 + q 12 √ −λ8 , b3 = − 1 2 ν ∓ θ 2 − q 4 √ −λ8 . (10) Substituting b1 = b = 0 and the values b2, b3 given by Eq. (10) into Eq. (5) leads to x = 1 3 √ −λ8  1 2 (ν ± θ) + q 4 √ −λ8   + 1 6 √ −λ8 υ˙ υ . (11) Note that the selection of a1 and d such that b1 = b = 0 provides the most useful algorithm. Subtracting side with side of Eqs. (4) and (5), respectively, we obtain (a2 − b2) x+ a3 − b3 = u˙ u − 1 2 υ˙ υ − x˙ x+ d . Substituting the value of x given by Eq. (11) into the above equation leads to 2 3 (ν ∓ θ)− 2 3 υ˙ υ = u˙ u − 1 2 υ˙ υ − x˙ x+ d . Integrating the equation with respect to t yields υ−2/3e2/3(ν∓θ)t = u (x+ d)υ1/2 . Taking a new function ξ into consideration as υ−2/3e2/3(ν∓θ)t = u (x+ d)υ1/2 = 1 ξ2 , (12) A coupling successive approximation method for solving Duffing equation . . . 81 therefore we obtain 1 2 υ˙ υ = 3 2 ξ˙ ξ + 1 2 (ν ∓ θ) . (13) Based on given Eq. (13) and Eq. (11), a transformation is obtained as x = − 2 3λ [ q + 4ν √ − λ 8 ] + 1 2 √ −λ 8 ξ˙ ξ . (14) 2.2. Transforming the initial equation to the resulting one Eq. (3) can be rewritten as d dt ( u˙ u + 1 2 υ˙ υ ) + u˙2 u2 − ( 1 2 υ˙ υ )2 = σ + p cosω t x+ d . (15) Multiplying both sides of the equation by u/υ1/2 leads to u υ1/2 d dt ( u˙ u + 1 2 υ˙ υ ) + u υ1/2 [( u˙ u )2 − ( 1 2 υ˙ υ )2] = (σ + p cosωt) u (x+ d)υ1/2 . Integrating the obtained equation with respect to t, we find an initial integral u υ1/2 ( u˙ u + 1 2 υ˙ υ ) = D2 + t∫ 0 (σ + p cosωt) u (x+ d)υ1/2 dt, (16) where D2 is an integral constant. Substituting the values a1 = 1, b = 0, parameters a2, a3 given by Eq. (10) and Eq. (14) into Eq. (4) leads to u˙ u = x˙ x+ d − 1 2 ξ˙ ξ + 1 2 (ν ∓ θ) (17) and integrating it with respect to t yields u = (x+ d) ξ−1/2e1/2(ν∓θ)t. (18) From Eqs. (12), (13), (17) and (18) one can get u υ1/2 ( u˙ u + 1 2 υ˙ υ ) = 1 ξ2 [ x˙+ (x+ d) ( ξ˙ ξ + ν ∓ θ )] . (19) Using Eqs. (8) and (14), the formulae of x˙ and x+ d are expressed as x˙ = 1 2 √ −λ8 ( ξ¨ ξ − ξ˙2 ξ2 ) , x + d = 1 2 √ −λ8 [ ξ˙ ξ + (ν ± θ) ] . 82 Dao Huy Bich, Nguyen Dang Bich Substituting the above-obtained formulae into Eq. (19) gives u υ1/2 ( u˙ u + 1 2 υ˙ υ ) = 1 2 √ −λ8 ξ 3 [ ξ¨ + 2νξ˙ − 1 2 Kξ ] , (20) where K = k − 4 3 q2 λ − 8 3 ν2. From Eqs. (12), (16), and (20), the resulting equation is obtained ξ¨ + 2νξ˙ − 1 2 Kξ = f (ξ, t) , (21) where f (ξ, t) = 2 √ − λ 8 ξ3  D2 + t∫ 0 (σ + p cosωt) 1 ξ2 dt   . The formula (14) and Eq. (21) are the transformation and the resulting equation that the present paper is looking for. In order to formulate a successive approximation method, the right hand side of Eq. (21) is written as f (ξ, t) = 2 √ −λ 8 η (ξ, t) ξ, (22) where η (ξ, t) = ξ2  D2 + t∫ 0 (σ + p cosω t) 1 ξ2 dt   . (23) 3. EQUATION SOLVED BY THE COUPLING SUCCESSIVE APPROXIMATION METHOD Finding an analytic approximated solution to Eq. (21) by the coupling successive approximation method is carried out by continuous loops of iteration. Each loop contains continuously iterative steps. 3.1. Loops of iteration In loop “0”th, we solve the linear differential equation (21) without the right hand side to find the solution of ξ0 (t). In the first loop, substituting ξ (t) = ξ0 (t) in the right hand side of Eq. (21) and solving the obtained linear differential equation we find ξ1 (t) and so on. In loop n−1th, the value ξn−1 (t) is found. The function η(ξn−1, t) is computed by the formula (22) η (ξn−1, t) = (ξn−1) 2  D2 + t∫ 0 (σ + p cosω t) 1 (ξn−1)2 dt   , (24) A coupling successive approximation method for solving Duffing equation . . . 83 the iteration scheme of successive approximation method is introduced as following ξ¨n + 2νξ˙n − 1 2 Kξn = 2 √ − λ 8 η(ξn−1, t)ξn−1, n = 1, 2, 3 . . . (25) By solving Eq. (25), where the right hand side is a known function and taking into account Eq. (9), the analytic approximated solution in nth loop of iteration is obtained ξn = √ −λ 8 θ yn−1 − √ −λ 8 θ zn−1 +D3e −(ν−θ)t −D4e −(ν+θ) t, (26) where yn−1 (t) = e −(ν−θ)t   t∫ 0 η(ξn−1, t)ξn−1e (ν−θ)tdt   , (27) zn−1 (t) = e −(ν+θ)t t∫ 0 η (ξn−1, t) ξn−1e (ν+θ)tdt, (28) θ = [ 1 2 ( k − 4 3 q2 λ − 2 3 ν2 )] 1 2 . Examining Eq. (26) we can predict some characteristics of solution. If k − 4 3 q2 λ − 2 3 ν2 > 0 then θ is a real number, the solution describes an oscillation depending only one excited frequency ω. If k − 4 3 q2 λ − 2 3 ν2 < 0 then θ is an imaginary number, i.e θ = iϕ with ϕ = [ 1 2 ( 2 3 ν2 + 4 3 q2 λ − k )]1/2 , where ϕ plays the role of a new frequency of a nonlinear vibration. The solution (26) describes a complex oscillation with many frequencies: excited frequency ω, vibration frequency ϕ and combined frequency of ω and ϕ, that the chaotic characteristics of solution may be predicted. Each function in the sequence ξ0(t), ξ1(t), . . . , ξn−1(t), ξn(t) can be determined from the one immediately proceeding it by solving the respective linear differential equation (25). The process is stopped when the condition max n ‖ξn (t)− ξn−1 (t)‖ < ε is achieved, where ε is a small positive number as required. But the convergence proving of this process is very complicated. Thus, a coupling successive approximation method based on Eq. (22) must be devel- oped with the iterative steps as following: in each loop of iteration, continuously iterative steps are carried out. 84 Dao Huy Bich, Nguyen Dang Bich 3.2. Iterative steps in each loop In the loop nth, when the iterative step mth is carried out, the value η (ξn−1, t) is fixed. This value is taken at the end of the previous loop (loop n− 1th). At this point, the iteration scheme of the coupling successive approximation method for the loop nth and the iterative step mth is expressed as ξ¨n,m + 2νξ˙n,m − 1 2 Kξn,m = 2 √ −λ 8 η (ξn−1, t) ξn,m−1, n = 1, 2, 3 . . . , m = 1, 2, 3 . . . (29) where n denotes the number of loop and n-the number of iterative step. The approximate solution ξn−1 (t) in the last loop n − 1 th is taken as an initial approximation at the iterative step “0”th of the loop nth, denoted as ξn,0 (t). Thus, that requires ξn−1 (t) = ξn,0 (t) . Solving Eq. (29), where the right hand side is a known function, we have ξn,m = √ −λ 8 θ yn,m−1 − √ −λ 8 θ zn,m−1 +D3e −(ν−θ) t −D4e −(ν+θ)t, (30) where yn,m−1 (t) = e −(ν−θ)t   t∫ 0 η(ξn,0, t)ξn,m−1e (ν−θ)tdt   , zn,m−1 (t) = e −(ν+θ)t t∫ 0 η (ξn,0, t) ξn,m−1e (ν+θ)tdt. (31) From which ξ˙n,m = − (ν − θ) √ −λ 8 θ yn,m−1 + (ν + θ) √ −λ 8 θ zn,m−1 − (ν − θ)D3e −(ν−θ) t + (ν + θ)D4e −(ν+θ) t, ξ¨n,m = (ν − θ) 2 √ −λ 8 θ yn,m−1 − (ν + θ) 2 √ −λ 8 θ zn,m−1 + (ν − θ) 2D3e −(ν−θ)t − (ν + θ)2D4e −(ν+θ)t + 2 √ −λ 8 η (ξn,0, t) ξn,m−1. (32) Remarks : If each loop is carried out with only one step, the coupling successive method will return to the single successive method as mentioned in Section 3.1. A coupling successive approximation method for solving Duffing equation . . . 85 4. FINDING INTEGRAL CONSTANT In the loop nth and the step m+ 1th, basing on Eq. (14), we have xn,m+1 = − 2 3λ ( 4ν √ −λ 8 + q ) + 1 2 √ −λ 8 ξ˙n,m+1 ξn,m+1 , (33) x˙n,m+1 = 1 2 √ −λ 8   ξ¨n,m+1 ξn,m+1 − ( ξ˙n,m+1 ξn,m+1 )2 . (34) Assumed that ξn,m+1 |t=0 = ξ 0 n,m+1, ξ˙n,m+1 |t=0 = ξ˙ 0 n,m+1, xn,m+1 |t=0 = x0, x˙n,m+1 |t=0 = x˙0. From the relations (30)-(32) ξ0n,m+1 = D3 −D4, ξ˙ 0 n,m+1 = − (ν − θ)D3 + (ν + θ)D4, ξ¨n,m+1 |t=0 = (ν − θ) 2D3 − (ν + θ) 2D4 + 2 √ −λ 8 η ( ξ0n,0, 0 ) ξ0n,m, according to Eq. (24) η(ξ0n,0, 0) = (ξ 0 n,0) 2D2. Thus, it can be written as ξ¨n,m+1 |t=0 = (ν − θ) 2D3 − (ν + θ) 2D4 + 2 √ −λ 8 ( ξon,0 )2 ξ0n,mD2. From the above relations x0 = − 2 3λ ( 4ν √ −λ 8 + q ) + 1 2 √ −λ 8 − (ν − θ)D3 + (ν + θ)D4 D3 −D4 , (35) x˙0 = 1 2 √ −λ8   (ν − θ)2D3 − (ν + θ) 2D4 + 2 √ −λ8 ( ξ0n,0 )2 ξ0n,mD2 D3 −D4 − [ − (ν − θ)D3 + (ν + θ)D4 D3 −D4 ]2} (36) Solving the Eqs. (35) and (36) in terms of the variables D2, D4 or D2, D3, here we find D2, D4 D4 = D3 x0 + 2 3λ ( 4ν √ −λ 8 + q ) + 1 2 q −λ 8 (ν − θ) x0 + 2 3λ ( 4ν √ −λ 8 + q ) + 1 2 q −λ 8 (ν + θ) , (37) 86 Dao Huy Bich, Nguyen Dang Bich D2 = 2θD3( ξ0n,0 )2 ξ0n,m { 2 √ −λ8 [ x0 + 2 3λ ( 4ν √ −λ 8 + q ) + ν + θ ]} .  x˙0 + 2 √ −λ 8 [ x0 + 2 3λ ( 4ν √ −λ 8 + q )]2 + (ν + θ)2 2θ  x0 + 2 3λ ( 4ν √ −λ 8 + q ) + 1 2 √ −λ 8 (ν − θ)   − (ν − θ)2 2θ  x0 + 2 3λ ( 4ν √ −λ 8 + q ) + 1 2 √ −λ 8 (ν + θ)     , (38) where D3 is an elective constant, it is possible to select D3 = 1. 5. CALCULATING STEPS Step 1: Providing input parameters k, q, ν, λ, p, ω, x0, x˙0. Step 2: Finding σ based on Eq. (7), θ based on Eq. (9) and other input parameters. Step 3: Finding the integral coefficients based on Eq. (37) and Eq. (38). Step 4: Getting the 0th approximation ξ1,0 = D3e −(ν−θ)t or ξ1,0 = D4e −(ν+θ)t, and deter- mining the initial condition ξ0n,0, ξ 0 n,m. Step 5: Finding ξn,m+1, ξ˙n,m+1 based on Eq. (30), Eq. (31) and Eq. (32). Step 6: Finding xn,m+1 (t) , x˙n,m+1 (t)based on Eq. (33) and Eq. (34). Step 7: Plotting {Re [xn,m+1 [t]] , Re [x˙n,m+1 [t]] , Im [xn,m+1 [t]]}. Step 8: Plotting {Re [xn,m+1 [t]] , Re [x˙n,m+1 [t]]}. Step 9: Plotting {Re [xn,m+1 [t]]}. The approximations in the first step of the first loop: ξ1,1, ξ˙1,1, ξ¨1,1 are found by integrating in closed form, next the approximations are found by integrating numerically. 6. APPLICATION AND ASSESSMENT OF SOLUTION PROPERTIES 6.1. Exact solution The resulting equation of the proposed method can be used to find exact solutions in some particular cases: Consider an initial equation x¨+ λx3 + kx = 0. It is derived from the Eq. (1) when ν = 0, q = 0, p = 0. In this case, the transformation (14) and the solving Eq. (21) can be written as x = 1 2 √ −λ 8 ξ˙ ξ , ξ¨ − 1 2 kξ = 0, A coupling successive approximation method for solving Duffing equation . . . 87 where: σ = 0 because ν = 0, q = 0 and selected constant can be chosen D2 = 0. Based on the resulting equations, we find ξ = C cosh( √ k 2 t+ φ), and x = √ −k λ tanh (√ k 2 t+ φ ) , (39) x˙ = k √ −1 2λ 1 cosh2 (√ k 2 t+ φ ) , (40) in which: C, φ - integral constants. -0.6 -0.4 -0.2 0.2 0.4 0.6 x t , mH L H L 0.1 0.2 0.3 0.4 0.5 ( ), ( / )x t m s& -10 -5 5 10 t, sH L 0.1 0.2 0.3 0.4 0.5 ( ), ( / )x t m s& -10 -5 5 10 t, sH L -0.6 -0.4 -0.2 0.2 0.4 0.6 x t , mH L H L Fig. 1. Comparison of exact solution (39) with numerical solutions using Mathematica 7 and the solution at first approximation of the proposed coupling successive approximation method (CSAM); k = 1, λ = −2, x[0] = 0, x˙[0] = 0.5 With the initial set of parameters: k = 1, λ = −2, φ = 0 and the initial condition x0 = x (t)|t=0 = 0, x˙0 = x˙ (t)|t=0 = 0.5, 88 Dao Huy Bich, Nguyen Dang Bich the exact solution (39), (40) can be found x = √ 1 2 tanh (√ 1 2 t ) , x˙ = 1 2 1 cosh2 (√ 1 2 t ) . Solving the initial equation with the set of parameters: k = 1, λ = −2, x0 = 0, x˙0 = 0.5 by the numerical method using Mathematica 7 and the coupling successive approximation method (CSAM) then comparing obtained results with the exact solutions (39), (40) we can see that they are coincided exactly (see Fig. 1). 6.2. Real-valued solution In this case, θ is real and λ is smaller than zero. In fact, consider Eq. (1) with given parameters k = 0.12, q = 0.0, λ= −1.0, ν = 0.4, ω = 1.0, p= 1.558, x0 = −0.4, x˙0 = −1.0. In this case θ can be evaluated as θ = 0.0816497. -0.10 -0.08 -0.06 -0.04 -0.02 x t , mH L H L -0.04 -0.02 0.02 0.04 0.06 ( ), ( / )x t m s& Fig. 2. Phase plane with t (50, 550), based on the results at the first approximation 600 700 800 900 1000 t, sH L -0.10 -0.08 -0.06 -0.04 -0.02 x t , mH L H L Fig. 3. Solution x(t), based on the results at the first approximation 600 700 800 900 1000 t, sH L -0.04 -0.02 0.02 0.04 0.06 ( ), ( / )x t m s& Fig. 4. Solution x˙(t), based on the results at the first approximation A coupling successive approximation method for solving Duffing equation . . . 89 Remarks: When θ is real and λ < 0, the obtained solution is real-valued one. The phase curves illustrated in Fig. 2 are smooth. They are not intersecting, inter- twined into a closed ring. The curve x(t) is periodic with decreasing amplitude, and the solution x(t) shows a stable motion (Figs. 3 and 4). 6.3. Complex-valued solutions Complex-valued solutions have two components, the real and imaginary, Re [x (t)] , Im [x (t)]. Differentiated complex-valued solution with respect to time also has two compo- nents, the real and imaginary, Re[x˙(t)], Im[x˙(t)]. From Eq. (1) and the equivalent Eq. (3), the initial integral (19) including the four components mentioned above is founded. There- fore only three components are independent. The three components form a phase space, which is different from a phase plane in the case of real valued solution [15]. Consider Eq. (1) with the following parameters k = −1/5, q = 0.0, λ = 8/15, ν = 1/50, ω = 0.32, p= 0.4, x0 = −0.4, x˙0 = −1.0. In this case θ can be evaluated as an imaginary number θ = 0.316493i, thus θ = iϕ where ϕ is real. Remarks : The coefficient of nonlinear term λ = 8/15 is twice as big as the coefficient of linear term k = −1/5. Thus, it is not suitable to use the assumption of small parameters in solving this problem. In this example, θ = 0.316493i, the solution is a complex-valued and chaotic one. From time period t(5000, 6000) to t(9500, 10500), the phase space has a relatively stable structure (Figs. 5 and 6). The curves in the phase space (Fig. 5) intersect. They are not smooth and have complex behavior. ()x té ùë û&Re (m/s) Im (m) ()x té ùë ûRe (m) x té ùë û() Fig. 5. Phase space with t(5000, 6000), based on the results at the first approximation ()x té ùë û&Re (m/s) Im (m) ()x té ùë ûRe (m) x té ùë û() Fig. 6. Phase space with t (9500,10500), based on the results at the first approximation The curves illustrating motion x(t) (Fig. 7) do not behave with any patterns, being affected by not only global motion, but also local motion. They never repeat themselves 90 Dao Huy Bich, Nguyen Dang Bich 5200 5400 5600 5800 6000 t, sH L -0.66 -0.64 -0.62 -0.60 -0.58 -0.56 Re x t , m@ H LD H L Fig. 7. The real component of solution x(t), based on the results at the first approximation Fig. 8. Poincaré section of the phase space in Fig. 5 with Im[x(t)] = 0 and there is no sign of resonance, although the frequency of exciting force ω = 0.32 is very close to the vibration frequency ϕ = 0.316493. Poincaré section (Fig. 8) includes a collection of points. Thus, the chaotic property of the solution in this example is proved. 6.4. Chaotic solution As can be seen that the indication of the chaotic solution to the Duffing equation is shown by the factor θ (see Eq. (9)), when θ = iϕ, ϕ is real number. Consider Eq. (1) with the following parameters k = 0.0, q = 0.0, λ = 1.0, ν = 0.025, ω = 1.0, p = 7.5, x0 = −0.4, x˙0 = −1.0. In this case θ can be evaluated as θ = 0.0144338i. Im (m) x té ùë û() ()x té ùë ûRe (m) ()x té ùë û&Re (m/s) Fig. 9. Phase space with t (150,1150), based on the results at the first approximation -0.10 -0.05 0.05 0.10 0.15 0.20 0.25 Re x t , m@ H LD H L -0.20 -0.15 -0.10 -0.05 0.05 0.10 0.15 Re x t , m@ H LD H L/s . Fig. 10. Phase plane with t (150,1150), based on the results at the first approximation A coupling successive approximation method for solving Duffing equation . . . 91 700 800 900 1000 1100 t, sH L -0.1 0.1 0.2 Re x t , m@ H LD H L Fig. 11. The real component of solution x(t), based on the results at the first approximation Fig. 12. Poincaré section of the phase space in Fig. 9 with Im[x(t)] = 0 Remarks: The coefficient λ = 1.0 and the coefficient of exciting force p = 7.5 have finite values different from zero and large, meanwhile the coefficient of the linear term k = 0. Thus, it is not suitable to use the assumption of small parameters in solving this problem. In this example, θ = 0.0144338i, the solution is complex-valued one. Based on the indication, θ = iϕ, the solution is chaotic. The curves in the phase space (Fig. 9), and the phase plane (Fig. 10) are rough, creased, intersecting and intertwined. The curves of the real-valued component of solution x(t) cluster together. They do not repeat each with other, but they have similar structure (Fig. 11). The clusters are thus considered sustainable. Poincaré section (Fig. 12) consists of a set of points. Thus, the chaotic property of the solution in this example is proved. 7. CONCLUSION Findings of the paper are summarized as follows: 1. An algorithm to solve the Duffing equation is proposed, in which a method to transform the initial equation to the resulting equation and a coupling successive approx- imation method (CSAM) to solve the resulting equation are presented. 2. Based on the proposed algorithm, the analytic approximated solutions obtained may be real-valued, complex-valued. 3. The indication of chaotic solutions to Duffing equation is found. The indication is k − 4 3 q2 λ − 2 3 ν2 < 0. 4. When the indication is satisfied, the Duffing equation has complex-valued solu- tions, and the phase curve is a spatial curve in a phase space in stead of the phase plane. The Poincare section consists of a set of points and the solution is chaotic. 92 Dao Huy Bich, Nguyen Dang Bich 5. A formula to compute the vibration frequency of Duffing equation is obtained ϕ = [ 1 2 ( 2 3 ν2 + 4 3 q2 λ − k )]1/2 , this frequency is only dependent on the parameters of the equation, and independent of parameters of exciting forces. 6. The structure of chaotic solutions is revealed, since the first approximation solu- tion is able to be expressed in an analytic form. The chaotic solution consists of complex vibration with many frequencies: exciting frequency ω, vibration frequency ϕ and com- bined frequency of ω and ϕ, with bounded amplitude. The chaotic solution is going towards the attracting set as the process towards limit. 7. Using the resulting equation of the proposed method an exact analytical solution can be found to some specific Duffing equations without right hand side. REFERENCES [1] A. Volmir. Non-linear dynamics of plates and shells. Science Edition M, (1972). [2] D. Evenson and R. Fulton. Some studies of the non-linear dynamic response of shell- type structures. In the book. Dynamic stability structures. Proc. of Int. confer. (ed. by Herrmann G.) Oxford, (1967), pp. 237–254. [3] F. Alijani, M. Amabili, K. Karagiozis, and F. Bakhtiari-Nejad. Nonlinear vibrations of functionally graded doubly curved shallow shells. Journal of Sound and Vibration, 330, (7), (2011), pp. 1432–1454. [4] B. Wu, C. Lim, and Y. Ma. Analytical approximation to large-amplitude oscillation of a non-linear conservative system. International Journal of Non-linear Mechanics, 38, (7), (2003), pp. 1037–1043. [5] D. H. Bich, D. V. Dung, and V. H. Nam. Nonlinear dynamical analysis of eccentrically stiffened functionally graded cylindrical panels. Composite Structures, 94, (8), (2012), pp. 2465–2473. [6] D. H. Bich and N. X. Nguyen. Nonlinear vibration of functionally graded circular cylindrical shells based on improved Donnell equations. Journal of Sound and Vibra- tion, 331, (25), (2012), pp. 5488–5501. [7] N. Dao, T. Chi, and N. Dung. Chaos phenomenon in a nonlinear Mathieu oscillator. In Proceeding of the Seventh National Congress on Mechanics, Hanoi, (December, 2002), pp. 18–20. [8] B. P. Flannery, W. H. Press, S. A. Teukolsky, and W. Vetterling. Numerical recipes in C. Press Syndicate of the University of Cambridge, New York, (1992). [9] A. K. Chopra. Dynamics of structures, Vol. 3. Prentice Hall New Jersey, (1995). [10] Wolfram.com/products/player/. [11] A. H. Nayfeh and B. Balachandran. Applied nonlinear dynamics: Analytical, compu- tational and experimental methods. John Wiley & Sons, (2008). [12] Y. A. Mitroposkii and N. V. Dao. Lectures on asymptotic methods on nonlinear dynamics. Vietnam National University Publishing House, (2003). REFERENCES 93 [13] P. Yu, Y. M. Desai, N. Popplewell, and A. H. Shah. The Krylov-Bogoliubov and Galerkin methods for nonlinear oscillations. Journal of Sound and Vibration, 192, (2), (1996), pp. 431–438. [14] G. Chen. Application of generalized Galerkin’s method to non-linear oscillations of two-degree-freedom systems. Journal of sound and vibration, 119, (2), (1997), pp. 225–242. [15] D. H. Bich and N. D. Bich. On the characteristics of solutions to the Duffing equations. In Proceeding of the IX National Congress on Mechanics, Hanoi, (December, 2012). pp. 58–66. [16] J. H. He, G. C. Wu, and F. Austin. The variational iteration method which should be followed. Nonlinear Science Letters A-Mathematics, Physics and Mechanics, 1, (1), (2010), pp. 1–30. [17] J. H. He and X. H. Wu. Construction of solitary solution and compaction-like solution by variational iteration method. Chaos Solitons & Fractals, 29, (2006). [18] J. H. He. Preliminary report on the energy balance for nonlinear oscillations. Me- chanics Research Communications, 29, (2), (2002), pp. 107–111. [19] T. O¨zis¸ and A. Yıldırım. Determination of the frequency-amplitude relation for a Duffing-harmonic oscillator by the energy balance method. Computers & Mathematics with Applications, 54, (7), (2007), pp. 1184–1187. [20] S.-Q. Wang and J.-H. He. Nonlinear oscillator with discontinuity by parameter- expansion method. Chaos, Solitons & Fractals, 35, (4), (2008), pp. 688–691. [21] S. Durmaz, S. A. Demirbag˘, and M. Kaya. High order Hamiltonian approach to nonlinear oscillators. International Journal of Nonlinear Sciences and Numerical Simulation, 11, (8), (2010), pp. 565–570. [22] J.-H. He. Hamiltonian approach to nonlinear oscillators. Physics Letters A, 374, (23), (2010), pp. 2312–2314. [23] J.-H. He. Variational approach for nonlinear oscillators. Chaos, Solitons & Fractals, 34, (5), (2007), pp. 1430–1439. [24] M. Kaya, S. Durmaz, and S. A. Demirbag. He’s variational approach to multiple cou- pled nonlinear oscillators. International Journal of Nonlinear Sciences and Numerical Simulation, 11, (10), (2010), pp. 859–866. [25] A. Kimiaeifar, A. Saidi, G. Bagheri, M. Rahimpour, and D. Domairry. Analytical solution for Van der Pol-Duffing oscillators. Chaos, Solitons & Fractals, 42, (5), (2009), pp. 2660–2666. VIETNAM ACADEMY OF SCIENCE AND TECHNOLOGY VIETNAM JOURNAL OF MECHANICS VOLUME 36, N. 2, 2014 CONTENTS Pages 1. Dao Huy Bich, Nguyen Dang Bich, A coupling successive approximation method for solving Duffing equation and its application. 77 2. Nguyen Thai Chung, Hoang Xuan Luong, Nguyen Thi Thanh Xuan, Dynamic stability analysis of laminated composite plate with piezoelectric layers. 95 3. Vu Le Huy, Shoji Kamiya, A direct evidence of fatigue damage growth inside silicon MEMS structures obtained with EBIC technique. 109 4. Nguyen Tien Khiem, Duong The Hung, Vu Thi An Ninh, Multiple crack identification in stepped beam by measurements of natural frequencies. 119 5. Nguyen Hong Son, Hoang Thi Bich Ngoc, Dinh Van Phong, Nguyen Manh Hung, Experiments and numerical calculation to determine aerodynamic char- acteristics of flows around 3D wings. 133 6. Gulshan Taj M. N. A., Anupam Chakrabarti, Mohammad Talha, Free vi- bration analysis of four parameter functionally graded plate accounting for realistic transverse shear mode. 145