On the analytical approximation of the quadratic non-linear oscillator by modified extended iteration method

Non-linear dynamic problems have fascinated applied mathematicians, physicists and engineers from a long time. Over the past few decades, applications in solid and structural mechanics as well as fluid mechanics have appeared, and currently, there is widespread interest in non-linear oscillators, strange attractors, as well as the chaotic and dynamical systems theories in the engineering and applied science communities.

Physical and mechanical oscillatory systems are often governed by non-linear differential equations. Unfortunately, with the exception of a few particular cases, the exact analytical solutions of such equations cannot be determined. In many cases, it is possible to replace the non-linear differential equation by a corresponding linear differential equation that approximates the original non-linear equation closely, to give useful results. Often such linearisation is not feasible or possible, and for this situation, the original differential equation itself must be directly dealt with.

However, in many cases, it is possible to compute accurate approximate analytical solutions of the equations. A large number of approximate methods, such as perturbation [1, 2, 3, 4, 5, 6], harmonic balance (HB) [7, 8, 9, 10, 11, 12, 13], homotopy perturbation [14], iteration [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28] and so on, are commonly used for solving non-linear oscillatory systems. The perturbation method is mainly used for small non-linear problems. On the other hand, HB and iteration methods are mostly used for strong and small non-linear problems.

One important class of non-linear oscillators is conservative oscillators, in which the restoring force is not dependent on time, the total energy is constant and any oscillation is stationary. Despite the great elegance and simplicity of such equations, the solutions of specific problems are significantly hard to derive. Finding innovative methods to analyse and solve these equations has become an interesting subject in the field of ordinary and partial differential equations and dynamical systems. Using non-linear equations for most real-life problems is not always possible, and sometimes, it is not even advantageous to express exact solutions of non-linear differential equations explicitly in terms of elementary functions or independent spatial and/or temporal variables; however, it is possible to find approximate solutions.

Perturbation means grossly small change; so, the method is adopted when the non-linearity is small. Thus, in case of strong non-linearities, the perturbation method is not generally adopted. It is used to construct a uniformly valid periodic solution to second-order non-linear differential equations. A critical feature of the technique is a middle step that breaks the problem into ‘solvable’ and ‘perturbation’ parts. Perturbation theory is applicable if the problem at hand cannot be solved exactly but can be formulated by adding a ‘small’ term to the mathematical description of the exactly solvable problem.

The HB method is a procedure for determining analytical approximations to the periodic solutions of differential equations by using a truncated Fourier series representation. An important advantage of the method is that it can be applied to non-linear oscillatory problems for which the non-linear terms are not ‘small’, i.e. no perturbation parameter need to exist. A disadvantage of the method is that it is difficult to predict, for a given non-linear differential equation, whether a first-order HB calculation will provide a sufficiently accurate approximation to the periodic solution.

The iterative technique is particularly useful for calculating approximate periodic solutions and the corresponding frequencies of truly non-linear oscillators for small and large amplitudes of oscillation.

The main intention of this research is to investigate the approximate analytical solutions using the modified extended iterative method to decompose the secular term, so that the solution can be obtained by the iterative procedure. This means that we can use the extended iterative method to investigate many non-linear problems. The main thrust of this technique is that the obtained solution rapidly converges to the exact solutions.

An extended iterative method is used to obtain the analytical solution of the quadratic non-linear oscillator. The procedure may be briefly described as follows.

A non-linear oscillator is modelled by the following expression: (1) x¨+f (x)=0, x (0)=A, x˙ (0)=0 , \ddot x + f{\kern 1pt} (x) = 0,{\kern 1pt} {\kern 1pt} x{\kern 1pt} (0) = A,{\kern 1pt} {\kern 1pt} \dot x{\kern 1pt} (0) = 0{\kern 1pt}, where the overdots denote differentiation with respect to time t.

First, we choose the natural frequency Ω of this system. Then, adding Ω²x on both sides of Eq. (1), we obtain (2) x¨+Ω2 x=Ω2 x−f(x)≡G(x, Ω). \ddot x + {\Omega ^{\rm{2}}}{\kern 1pt} x = {\Omega ^{\rm{2}}}{\kern 1pt} x - f(x) \equiv G(x,{\kern 1pt} \Omega). The extended iterative scheme is (3) x¨k+1+Ωk2 xk+1=G(xk−1,Ω)+Gx(xk−1,Ω)(xk−xk−1); k=1,2,... {\ddot x_{k + 1}} + \Omega _k^2{\kern 1pt} {x_{k + 1}} = G({x_{k - 1}},\Omega) + {G_x}({x_{k - 1}},\Omega)({x_k} - {x_{k - 1}});{\kern 1pt} {\kern 1pt} {\kern 1pt} k = 1,2,... where Gx=∂ G∂ x {G_x} = {{\partial {\kern 1pt} G} \over {\partial {\kern 1pt} x}} .

The right-hand side of Eq. (3) is essentially the first term in a Taylor series expansion of the function G(xk,x˙k) G({x_k},{\dot x_k}) at the point (xk−1,x˙k−1) ({x_{k - 1}},{\dot x_{k - 1}}) [29].

We have the direct iteration scheme of Eq. (2), as shown in Eq. (4): (4) x¨k+1+Ωk2xk+1=G(xk,Ωk); k=0, 1, 2, ... {\ddot x_{k + 1}} + \Omega _k^2{x_{k + 1}} = G({x_k},{\Omega _k});{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} k = 0,{\kern 1pt} 1,{\kern 1pt} 2,{\kern 1pt}... and x_k+1 satisfies the condition (5) xk+1(0)=A. {x_{k + 1}}(0) = A. The initial estimate is considered to be the following [10]: (6) x0(t)=A cos (Ω0 t). {x_0}(t) = A{\kern 1pt} \cos {\kern 1pt} ({\Omega _0}{\kern 1pt} t). The above procedure gives the sequence of following solutions: x₁(t), x₂(t), x₃(t),.... The method can be extrapolated to any order of approximation; but due to growing algebraic complexity, the solution is confined to a lower order, usually the second [15].

Let us consider the following non-linear inverse oscillator: (7) x¨+x2 = 0 \ddot x + {x^2}{\kern 1pt} = {\kern 1pt} {\kern 1pt} 0 Adding Ω²x on both sides of Eq. (7), we get (8) x¨+Ω2x=Ω2x−x2=G(x,Ω) \ddot x + {\Omega ^{\rm{2}}}x = {\Omega ^{\rm{2}}}x - {x^2} = G(x,\Omega) where G(x,Ω) = Ω²x − x² and G_x(x,Ω) = Ω² − 2x.

According to Eq. (4), the direct iterative scheme of Eq. (8) is (9) x¨k+1+Ωk2xk+1=Ωk2xk−xk2 {\ddot x_{k + 1}} + \Omega _k^2{x_{k + 1}} = \Omega _k^2{x_k} - x_k^2 The first approximation x₁(t) and the frequency Ω₀ are obtained by substituting k = 0 in Eq. (9), and using Eq. (6), we get (10) x¨1+Ω02 x1=Ω02 x0−x02 {\ddot x_1} + \Omega _0^2{\kern 1pt} {x_1} = \Omega _0^2{\kern 1pt} {x_0} - x_0^2{\kern 1pt} where x₀(t) = A cos(Ω₀t) = A cos θ and θ = Ω₀t.

Now, substituting and expanding the right-hand side in a Fourier cosine series, Eq. (10) reduces to (11) x¨1+Ω02 x1=Ω02 A cosθ−(A cosθ)2x¨1+Ω02 x1=Ω02 A cosθ−A2 (0.848826 cos θ+0.169765 cos3θ −0.024252 cos5θ+0.008084 cos7θ−0.003675 cos9θ +0.001979 cos11θ)=(Ω02 −0.848826A)A cosθ−A2(0.169765 cos3θ+0.024252 cos5θ −0.008084 cos7θ+0.003675 cos9θ−0.001979 cos11θ ). \matrix{{{{\ddot x}_1} + \Omega _0^2{\kern 1pt} {x_1}} \hfill & {= \Omega _0^2{\kern 1pt} A{\kern 1pt} \cos \theta - {{(A{\kern 1pt} \cos \theta)}^2}} \hfill \cr {{{\ddot x}_1} + \Omega _0^2{\kern 1pt} {x_1}} \hfill & {= \Omega _0^2{\kern 1pt} A{\kern 1pt} \cos \theta - {A^2}{\kern 1pt} (0.848826{\kern 1pt} \cos {\kern 1pt} \theta + 0.169765{\kern 1pt} \cos 3\theta} \hfill \cr {} \hfill & {\,\,\,\,\,\, - 0.024252{\kern 1pt} \cos 5\theta + 0.008084{\kern 1pt} \cos 7\theta - 0.003675{\kern 1pt} \cos 9\theta} \hfill \cr {} \hfill & {\,\,\,\,\,\, + 0.001979{\kern 1pt} \cos 11\theta)} \hfill \cr {} \hfill & {= (\Omega _0^2{\kern 1pt} - 0.848826A)A{\kern 1pt} \cos \theta - {A^2}(0.169765{\kern 1pt} \cos 3\theta + 0.024252{\kern 1pt} \cos 5\theta} \hfill \cr {} \hfill & {\,\,\,\,\, - 0.008084{\kern 1pt} \cos 7\theta + 0.003675{\kern 1pt} \cos 9\theta - 0.001979{\kern 1pt} \cos 11\theta {\kern 1pt}).} \hfill \cr} To avoid secular terms in the solution, we have to remove cos θ from the right-hand side of Eq. (11). Thus, we have (12) Ω02 A−0.848826A2=0 , Ω02=0.848826A2A, Ω0=0.921318A. \Omega _0^2{\kern 1pt} A - 0.848826{A^2} = 0{\kern 1pt},{\kern 1pt} {\kern 1pt} \Omega _0^2 = {{0.848826{A^2}} \over A},{\kern 1pt} {\kern 1pt} {\Omega _0} = 0.921318\sqrt A. This is the first approximate frequency of the oscillator. Note that Ωexact(A)=0.914681A {\Omega _{{\rm{exact}}}}(A) = 0.914681\sqrt A . After simplification, Eq. (11) reduces to (13) x¨1+Ω02 x1=A2(−0.169765 cos3θ+0.024252 cos5θ−0.008084 cos7θ +0.003675 cos9θ−0.001979 cos11θ). \matrix{{{{\ddot x}_1} + \Omega _0^2{\kern 1pt} {x_1} = {A^2}(- 0.169765\cos 3\theta + 0.024252{\kern 1pt} \cos 5\theta - 0.008084{\kern 1pt} \cos 7\theta} \hfill \cr {{\kern 70pt} + 0.003675{\kern 1pt} \cos 9\theta - 0.001979{\kern 1pt} \cos 11\theta).} \hfill \cr} The particular solution, x₁^(p)(t), is as follows: (14) x1(p)(t)=−0.169765A2−9Ω02+Ω02 cos3θ+0.024252A2−25Ω02+Ω02 cos5 θ+−0.008084A2−49Ω02+Ω02 cos7θ +0.003675A2−81Ω02+Ω02 cos9θ+0.001979A2−121Ω02+Ω02 cos11θ=A2Ω02(0.1697658 cos3θ−0.02425224 cos5 θ+0.00808448 cos7θ −0.00367580 cos9θ−0.001979120 cos11θ)=A(0.025cos3θ−0.001190cos5θ+0.000199cos7θ −0.000054cos9θ−0.000019cos11θ) \matrix{{x_1^{(p)}(t)} \hfill & {= {{- 0.169765{A^2}} \over {- 9\Omega _0^2 + \Omega _0^2}}{\kern 1pt} \cos 3\theta + {{0.024252{A^2}} \over {- 25\Omega _0^2 + \Omega _0^2}}{\kern 1pt} \cos 5{\kern 1pt} \theta + {{- 0.008084{A^2}} \over {- 49\Omega _0^2 + \Omega _0^2}}{\kern 1pt} \cos 7\theta} \hfill \cr {} \hfill & {\,\,\,\,\,\, + {{0.003675{A^2}} \over {- 81\Omega _0^2 + \Omega _0^2}}{\kern 1pt} \cos 9\theta + {{0.001979{A^2}} \over {- 121\Omega _0^2 + \Omega _0^2}}{\kern 1pt} \cos 11\theta} \hfill \cr {} \hfill & {= {{{A^2}} \over {\Omega _0^2}}\left({{{0.169765} \over 8}{\kern 1pt} \cos 3\theta - {{0.024252} \over {24}}{\kern 1pt} \cos 5{\kern 1pt} \theta + {{0.008084} \over {48}}{\kern 1pt} \cos 7\theta} \right.} \hfill \cr {} \hfill & {\,\,\,\left. {\,\,\, - {{0.003675} \over {80}}{\kern 1pt} \cos 9\theta - {{0.001979} \over {120}}{\kern 1pt} \cos 11\theta} \right)} \hfill \cr {} \hfill & {= A(0.025\cos 3\theta - 0.001190\cos 5\theta + 0.000199\cos 7\theta} \hfill \cr {} \hfill & {\,\,\, - 0.000054\cos 9\theta - 0.000019\cos 11\theta)} \hfill \cr} Therefore, the complete solution is (15) x1(t)=B1 cosθ+0.025Acos3θ−0.001190Acos5θ+0.000199Acos7θ −0.000054Acos9θ−0.000019Acos11θ . \matrix{{{x_1}(t)} \hfill & {= {B_1}{\kern 1pt} \cos \theta + {\rm{0}}{\rm{.025}}A{\rm{cos3}}\theta - {\rm{0}}{\rm{.001190}}A{\rm{cos5}}\theta + {\rm{0}}{\rm{.000199}}A{\rm{cos7}}\theta} \hfill \cr {} \hfill & {\,\,\,\,\, - 0.000054A{\rm{cos9}}\theta - {\rm{0}}{\rm{.000019}}A{\rm{cos11}}\theta {\kern 1pt}.} \hfill \cr} Using x₁(0) = A, we have B₁ = 0.976066A. Then, we obtain (16) x1(t)=A(0.976066 cosθ+0.025 cos3θ−0.00119048 cos5θ+0.00198413 cos7θ −0.0000541126 cos9θ−0.000019425 cos11θ). \eqalign{& {x_1}(t) = A(0.976066\cos \theta + 0.025\cos 3\theta - 0.00119048\cos 5\theta + 0.00198413\cos 7\theta \cr & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, - 0.0000541126\cos 9\theta - 0.000019425\cos 11\theta). \cr} This is the first approximate solution of the oscillator.

According to Eq. (3), the extended iterative scheme of Eq. (8) is (17) x¨k+1+Ωk2 xk+1= (Ω2 xk−xk2)+(Ω2−2xk) (xk−xk−1). =xk2+Ω2xk−1−2xkxk−1 \matrix{{{\kern 1pt} {{\ddot x}_{k + 1}} + \Omega _k^2{\kern 1pt} {x_{k + 1}}} \hfill & {= ({\Omega ^2}{\kern 1pt} {x_k} - x_k^2) + ({\Omega ^{\rm{2}}} - 2{x_k}){\kern 1pt} ({x_k} - {x_{k - 1}}).} \hfill \cr {\kern 1pt} \hfill & {= x_k^2 + {\Omega ^2}{x_{k - 1}} - 2{x_k}{x_{k - 1}}} \hfill \cr} The second approximation x₂(t) and the frequency Ω₁ are obtained by substituting k = 1 in Eq. (17), and using Eq. (7), we get Eq. (18): (18) x¨2+Ω12 x2=x12 +Ω2x0−2x1x0 {\ddot x_2} + \Omega _1^2{\kern 1pt} {x_2} = x_1^2{\kern 1pt} + {\Omega ^2}{x_0} - 2{x_1}{x_0} where x₀(t) and x₁(t) are given by Eqs (6) and (16).

Now substituting x₀(t) and x₁(t) and expanding the right-hand side in a Fourier cosine series, Eq. (17) reduces to (19) x¨2+Ω12 x2=(Ω12 0.976066A−0.816744A2)cosθ+(Ω12 0.025000A −0.193883A2)cos3θ+(−Ω12 0.001190A+0.014431A2)cos5θ +(Ω12 0.000198A−0.005635A2)cos7θ+(−Ω120.000054A +0.002713A2)cos9θ+(−Ω12 0.000019A−0.001441A2)cos11θ +0.000911A2cos13θ \matrix{{{{\ddot x}_2} + \Omega _1^2{\kern 1pt} {x_2} = (\Omega _1^2{\kern 1pt} 0.976066A - 0.816744{A^2})\cos \theta + (\Omega _1^2{\kern 1pt} 0.025000A{\kern 1pt}} \hfill \cr {{\kern 65pt} - 0.193883{A^2})\cos 3\theta + (- \Omega _1^2{\kern 1pt} 0.001190A + 0.014431{A^2})\cos 5\theta} \hfill \cr {{\kern 65pt} + (\Omega _1^2{\kern 1pt} 0.000198A - 0.005635{A^2})\cos 7\theta + (- \Omega _1^20.000054A} \hfill \cr {{\kern 65pt} + 0.002713{A^2})\cos 9\theta + (- \Omega _1^2{\kern 1pt} 0.000019A - 0.001441{A^2})\cos 11\theta} \hfill \cr {{\kern 65pt} + 0.000911{A^2}\cos 13\theta} \hfill \cr} To avoid secular terms in the solution, we have to remove cos θ from the right-hand side of Eq. (19). Thus, we have (20) Ω120.976066A−0.816744A2=0, Ω1=0.914752A. \Omega _1^20.976066A - 0.816744{A^2} = 0,\;\;{\Omega _1} = 0.914752\sqrt A. The particular solution, x₂^(p)(t), is (21) x2(p)(t)=1−9 Ω12+Ω12(Ω12 0.025000A−0.193883A2)cos3θ+ 1−25 Ω12+Ω12 (−Ω120.001190A+0.014431A2) cos5θ+1−48 Ω12+Ω12(Ω120.000198A −0.005635A2)cos7θ+1−81 Ω12+Ω12(−Ω120.000054A+0.002712A2) cos9θ+1−121 Ω12+Ω12(−Ω120.000019A −0.001441A2)cos11θ +1−169Ω12+Ω12(0.000911A2)cos13θ =−0.025000A8 cos3θ +0.193883A28Ω12cos3θ− 0.001190A24 cos5θ −0.014431A224Ω12 cos5θ−0.000198A48 cos7θ+0.005635A248 Ω12cos7θ +0.000054A80cos9θ−0.002712A280 Ω12cos9θ+0.000019A120cos11θ +0.001441A2120 Ω12cos11θ −0.000911A2168Ω12cos13θ =0.025838Acos 3θ−0.000669A cos5θ+0.000136A cos7θ −0.000040A cos9θ+0.000015A cos11θ −6.48356 ×10−6 A cos13θ \matrix{{x_2^{(p)}(t) = {1 \over {- 9{\kern 1pt} \Omega _1^2 + \Omega _1^2}}(\Omega _1^2{\kern 1pt} 0.025000A - 0.193883{A^2})\cos 3\theta + {\kern 1pt} {1 \over {- 25{\kern 1pt} \Omega _1^2 + \Omega _1^2}}{\kern 1pt}} \hfill \cr {{\kern 25pt}(- \Omega _1^20.001190A + 0.014431{A^2}){\kern 1pt} \cos 5\theta + {1 \over {- 48{\kern 1pt} \Omega _1^2 + \Omega _1^2}}(\Omega _1^20.000198A} \hfill \cr {{\kern 25pt} - 0.005635{A^2})\cos 7\theta + {1 \over {- 81{\kern 1pt} \Omega _1^2 + \Omega _1^2}}(- \Omega _1^20.000054A + 0.002712{A^2})} \hfill \cr {{\kern 25pt} \cos 9\theta + {1 \over {- 121{\kern 1pt} \Omega _1^2 + \Omega _1^2}}(- \Omega _1^20.000019A{\kern 1pt} - 0.001441{A^2})\cos 11\theta {\kern 1pt}} \hfill \cr {{\kern 20pt} + {1 \over {- 169\Omega _1^2 + \Omega _1^2}}(0.000911{A^2})\cos 13\theta} \hfill \cr {{\kern 20pt} = - {{0.025000A} \over {8{\kern 1pt}}}\cos 3\theta {\kern 1pt} + {{0.193883{A^2}} \over {8\Omega _1^2}}\cos 3\theta - {\kern 1pt} {{0.001190A} \over {24{\kern 1pt}}}\cos 5\theta {\kern 1pt}} \hfill \cr {{\kern 20pt} - {{0.014431{A^2}} \over {24\Omega _1^2}}{\kern 1pt} \cos 5\theta - {{0.000198A} \over {48{\kern 1pt}}}\cos 7\theta + {{0.005635{A^2}} \over {48{\kern 1pt} \Omega _1^2}}\cos 7\theta} \hfill \cr {{\kern 20pt} + {{0.000054A} \over {80}}\cos 9\theta - {{0.002712{A^2}} \over {80{\kern 1pt} \Omega _1^2}}\cos 9\theta + {{0.000019A} \over {120}}\cos 11\theta {\kern 1pt}} \hfill \cr {\kern 20pt} \matrix{+ {{0.001441{A^2}} \over {120{\kern 1pt} \Omega _1^2}}\cos 11\theta {\kern 1pt} - {{0.000911{A^2}} \over {168\Omega _1^2}}\cos 13\theta \hfill \cr \matrix{{= 0.025838A{\rm{cos}}{\kern 1pt} 3\theta - 0.000669A{\kern 1pt} {\rm{cos}}5\theta + 0.000136A{\kern 1pt} {\rm{cos}}7\theta} \hfill \cr {- 0.000040A{\kern 1pt} {\rm{cos}}9\theta + 0.000015A{\kern 1pt} {\rm{cos}}11\theta} \hfill \cr {- 6.48356{\kern 1pt} \times {{10}^ -}6{\kern 1pt} A{\kern 1pt} {\rm{cos}}13\theta {\kern 1pt}} \hfill \cr} \hfill \cr} \hfill \cr} Therefore, the complete solution is (22) x2(t)=B2 cosθ+0.025838Acos3θ−0.000669Acos5θ+0.000136Acos7θ −0.000040Acos9θ+0.000015Acos11θ−6.48356×10−6Acos13θ \matrix{{{x_2}(t) = {B_2}{\kern 1pt} {\rm{cos}}\theta + 0.025838A{\rm{cos}}3\theta - 0.000669A{\rm{cos}}5\theta + 0.000136A{\rm{cos}}7\theta} \hfill \cr {{\kern 40pt} - 0.000040A{\rm{cos}}9\theta + 0.000015A{\rm{cos}}11\theta - 6.48356 \times {{10}^{- 6}}A{\rm{cos}}13\theta} \hfill \cr} Using x₂(0) = A, we have B₂ = 0.974727A. Then, we obtain (23) x2(t)=0.974727Acosθ+0.0258379Acos3θ−0.000669Acos5θ +0.000136Acos7θ−0.000040Acos9 θ+0.000015Acos11θ −6.48356×10−6Acos13θ \matrix{{{x_2}(t) = 0.974727A{\rm{cos}}\theta + 0.0258379A{\rm{cos}}3\theta - 0.000669A{\rm{cos}}5\theta} \hfill \cr {{\kern 40pt} + 0.000136A{\rm{cos}}7\theta - 0.000040A{\rm{cos}}9{\kern 1pt} \theta + 0.000015A{\rm{cos}}11\theta} \hfill \cr {{\kern 40pt}- 6.48356 \times {{10}^{- 6}}A{\rm{cos}}13\theta} \hfill \cr} This is the second approximate solution of the oscillator.

Proceeding to the third level of iteration, satisfies the following equation: (24) x¨3+Ω22 x3=x22+Ω2x1−2x2x1 {\ddot x_3} + \Omega _2^2{\kern 1pt} {x_3} = x_2^2 + {\Omega ^2}{x_1} - 2{x_2}{x_1} where x1(t)=A(0.976066 cosθ+0.025 cos3θ−0.00119048 cos5θ +0.000198413 cos7θ−0.0000541126 cos9θ−0.000019425 cos11θ) \eqalign{& {x_1}(t) = A(0.976066\cos \theta + 0.025\cos 3\theta - 0.00119048\cos 5\theta \cr & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, + 0.000198413\cos 7\theta - 0.0000541126\cos 9\theta - 0.000019425\cos 11\theta) \cr} and x2(t)=A(0.974727 cosθ+0.0258379 cos3θ−0.000669 cos5θ+0.000136 cos7θ −0.000040 cos9θ+0.000015 cos11θ−6.48356×10−6cos13θ) \matrix{{{x_2}(t) = A(0.974727{\rm{cos}}\theta + 0.0258379{\rm{cos}}3\theta - 0.000669{\rm{cos}}5\theta + 0.000136{\rm{cos}}7\theta} \hfill \cr {{\kern 40pt} - 0.000040{\rm{cos}}9\theta + 0.000015{\rm{cos}}11\theta - 6.48356 \times {{10}^{- 6}}{\rm{cos}}13\theta)} \hfill \cr} Now substituting x₁(t) and x₂(t) and expanding the right-hand side in a Fourier cosine series, Eq. (24) reduces to (25) x¨3+Ω22x3=(−Ω220.974727A−0.815477A2)cosθ+(Ω220.025838A−0.193801A2)cos3θ+(−Ω220.000669A+0.013352A2)cos5θ+(Ω220.000136A−0.0058897A2)cos7θ+(−Ω220.0000399A+0.002779A2)cos9θ+(Ω220.000015A−0.001523A2)cos11θ+(−Ω220.000006A+0.000924A2)cos13θ−0.00597A2cos15θ+0.00412A2cos17θ−0.000296A2cos19θ+0.000219A2cos21θ+0.000167A2cos23θ+0.0001298A2cos25θ−0.000103A2cos27θ \matrix{{{{\ddot x}_3} + \Omega _2^2{x_3} =} \hfill & {(- \Omega _2^20.974727A - 0.815477{A^2})\cos \theta + (\Omega _2^20.025838A} \hfill \cr {} \hfill & {- 0.193801{A^2})\cos 3\theta + (- \Omega _2^20.000669A + 0.013352{A^2})\cos 5\theta} \hfill \cr {} \hfill & {+ (\Omega _2^20.000136A - 0.0058897{A^2})\cos 7\theta + (- \Omega _2^20.0000399A} \hfill \cr {} \hfill & {+ 0.002779{A^2})\cos 9\theta + (\Omega _2^20.000015A - 0.001523{A^2})\cos 11\theta} \hfill \cr {} \hfill & {+ (- \Omega _2^20.000006A + 0.000924{A^2})\cos 13\theta - 0.00597{A^2}\cos 15\theta} \hfill \cr {} \hfill & {+ 0.00412{A^2}\cos 17\theta - 0.000296{A^2}\cos 19\theta + 0.000219{A^2}\cos 21\theta} \hfill \cr {} \hfill & {+ 0.000167{A^2}\cos 23\theta + 0.0001298{A^2}\cos 25\theta - 0.000103{A^2}\cos 27\theta} \hfill \cr} To avoid secular terms in the solution, we have to remove from the right-hand side of Eq. (25). Thus, we have (26) −Ω220.974727A−0.815477A2=0, Ω2=0.91467A. - \Omega _2^20.974727A - 0.815477{A^2} = 0,\;\;{\Omega _2} = 0.91467\sqrt A. The particular solution, x3(p)(t) x_3^{(p)}(t) , is (27) x3(p)(t)=−18 Ω22(Ω22 0.025838A−0.193801A2)cos3θ − 124 Ω22 (−Ω220.000669A+0.013352A2) cos5θ−148 Ω22(Ω220.000136A −0.0058998A2)cos7θ−180 Ω22(−Ω220.0000399A+0.002779A2) cos9θ −1120 Ω22(Ω220.000015A−0.001523A2)cos11θ−1168Ω22 (−Ω22 0.000006A+0.000924A2)cos13θ−1224Ω22(−0.000597A2) cos15θ−1288 Ω22(0.000412A2)cos17θ−1360 Ω22(−0.000296A2) cos19θ −1440 Ω22(0.000219A2)cos21θ−1528 Ω22 (0.000167A2) cos23θ−1624 Ω22(0.0001298A2)cos25θ−1728 Ω22(−0.000103A2)cos27θ =0.025726Acosθ−0.000637Acos5θ+0.00014Acos7θ −0.000041Acos9θ+0.000015Acos11θ−6.53694×10−6Acos13θ +3.18798×10−6Acos15θ−1.71121×10−6Acos17θ +9.381306×10−7Acos19θ−5.94814×10−7Acos21θ −3.77331×10−7Acos23θ−2.4863×10−7Acos25θ +1.69177×10−7Acos27θ \matrix{{x_3^{(p)}(t) = - {1 \over {8{\kern 1pt} \Omega _2^2}}(\Omega _2^2{\kern 1pt} 0.025838A - 0.193801{A^2})\cos 3\theta {\kern 1pt} - {\kern 1pt} {1 \over {24{\kern 1pt} \Omega _2^2}}{\kern 1pt} {\kern 1pt}} \hfill \cr {{\kern 40pt} (- \Omega _2^20.000669A + 0.013352{A^2}){\kern 1pt} \cos 5\theta - {1 \over {48{\kern 1pt} \Omega _2^2}}(\Omega _2^20.000136A} \hfill \cr {{\kern 40pt} - 0.0058998{A^2})\cos 7\theta - {1 \over {80{\kern 1pt} \Omega _2^2}}(- \Omega _2^20.0000399A + 0.002779{A^2})} \hfill \cr {{\kern 40pt} \cos 9\theta {\kern 1pt} - {1 \over {120{\kern 1pt} \Omega _2^2}}(\Omega _2^20.000015A - 0.001523{A^2})\cos 11\theta - {1 \over {168\Omega _2^2}}} \hfill \cr {{\kern 40pt} (- \Omega _2^2{\kern 1pt} 0.000006A + 0.000924{A^2})\cos 13\theta - {1 \over {224\Omega _2^2}}(- 0.000597{A^2})} \hfill \cr {{\kern 40pt} \cos 15\theta - {1 \over {288{\kern 1pt} \Omega _2^2}}(0.000412{A^2})\cos 17\theta - {1 \over {360{\kern 1pt} \Omega _2^2}}(- 0.000296{A^2})} \hfill \cr {{\kern 40pt} \cos 19\theta {\kern 1pt} - {1 \over {440{\kern 1pt} \Omega _2^2}}(0.000219{A^2})\cos 21\theta - {1 \over {528{\kern 1pt} \Omega _2^2}}{\kern 1pt} (0.000167{A^2})} \hfill \cr \matrix{{\kern 40pt} \cos 23\theta - {1 \over {624{\kern 1pt} \Omega _2^2}}(0.0001298{A^2})\cos 25\theta - {1 \over {728{\kern 1pt} \Omega _2^2}}(- 0.000103{A^2})\cos 27\theta \hfill \cr {\kern 40pt} = 0.025726A\cos \theta - 0.000637A\cos 5\theta + 0.00014A\cos 7\theta \hfill \cr {\kern 45pt} - 0.000041A\cos 9\theta + 0.000015A\cos 11\theta - 6.53694 \times {10^{- 6}}A\cos 13\theta \hfill \cr {\kern 45pt} + 3.18798 \times {10^{- 6}}A\cos 15\theta - 1.71121 \times {10^{- 6}}A\cos 17\theta \hfill \cr {\kern 45pt} + 9.381306 \times {10^{- 7}}A\cos 19\theta - 5.94814 \times {10^{- 7}}A\cos 21\theta \hfill \cr {\kern 45pt} - 3.77331 \times {10^{- 7}}A\cos 23\theta - 2.4863 \times {10^{- 7}}A\cos 25\theta \hfill \cr {\kern 45pt} + 1.69177 \times {10^{- 7}}A\cos 27\theta \hfill \cr} \hfill \cr} Therefore, the complete solution is (28) x3(t)=B3cosθ+0.0257261Acos3θ−0.000637Acos5θ+0.000144Acos7θ−0.0000410252Acos9θ+0.0000150452Acos11θ−6.53694×10−6Acos13θ+3.18798×10−6Acos15θ−1.711121×10−6Acos17θ+9.81306×10−7Acos19θ−5.94814×10−7Acos21θ−3.77331×10−7Acos23θ−2.4863×10−7Acos25θ+1.69177×10−7Acos27θ \matrix{{{x_3}(t)} \hfill & {= {B_3}\cos \theta + 0.0257261A\cos 3\theta - 0.000637A\cos 5\theta} \hfill \cr {} \hfill & {+ 0.000144A\cos 7\theta - 0.0000410252A\cos 9\theta} \hfill \cr {} \hfill & {+ 0.0000150452A\cos 11\theta - 6.53694 \times {{10}^{- 6}}A\cos 13\theta} \hfill \cr {} \hfill & {+ 3.18798 \times {{10}^{- 6}}A\cos 15\theta - 1.711121 \times {{10}^{- 6}}A\cos 17\theta} \hfill \cr {} \hfill & {+ 9.81306 \times {{10}^{- 7}}A\cos 19\theta - 5.94814 \times {{10}^{- 7}}A\cos 21\theta} \hfill \cr {} \hfill & {- 3.77331 \times {{10}^{- 7}}A\cos 23\theta - 2.4863 \times {{10}^{- 7}}A\cos 25\theta} \hfill \cr {} \hfill & {+ 1.69177 \times {{10}^{- 7}}A\cos 27\theta} \hfill \cr} Using x₃(0) = A, we have B₃ = 0.974798A. Then, we obtain (29) x3(t)=0.974798Acosθ+0.025726Acos3θ−0.000637Acos5θ+0.000144Acos7θ−0.000041Acos9θ+0.000015Acos11θ−6.53694×10−6Acos13θ+3.18798×10−6Acos15θ−1.711121×10−6Acos17θ+9.81306×10−7Acos19θ−5.94814×10−7Acos21θ−3.77331×10−7Acos23θ−2.4863×10−7Acos25θ+1.69177×10−7Acos27θ \matrix{{{x_3}(t)} \hfill & {= 0.974798A\cos \theta + 0.025726A\cos 3\theta - 0.000637A\cos 5\theta} \hfill \cr {} \hfill & {+ 0.000144A\cos 7\theta - 0.000041A\cos 9\theta} \hfill \cr {} \hfill & {+ 0.000015A\cos 11\theta - 6.53694 \times {{10}^{- 6}}A\cos 13\theta} \hfill \cr {} \hfill & {+ 3.18798 \times {{10}^{- 6}}A\cos 15\theta - 1.711121 \times {{10}^{- 6}}A\cos 17\theta} \hfill \cr {} \hfill & {+ 9.81306 \times {{10}^{- 7}}A\cos 19\theta - 5.94814 \times {{10}^{- 7}}A\cos 21\theta} \hfill \cr {} \hfill & {- 3.77331 \times {{10}^{- 7}}A\cos 23\theta - 2.4863 \times {{10}^{- 7}}A\cos 25\theta} \hfill \cr {} \hfill & {+ 1.69177 \times {{10}^{- 7}}A\cos 27\theta} \hfill \cr} This is the third approximate solution of the oscillator.

Thus, Ω₀ , Ω₁ , Ω₂ , respectively obtained from Eqs (12), (20) and (26), represent the approximation of the frequencies, and x₁(t), x₂(t), x₃(t), respectively obtained from Eqs (16), (23) and (29), represent the corresponding approximate solutions of the oscillator represented by Eq. (7).

Let us consider the non-linear oscillator (30) x¨+x=−xx˙2 \ddot x + x = - x{\dot x^2} with initial condition x(x) = A, x˙(0)=0 \dot x(0) = 0 .

Obviously, Eq. (30) can be written as (31) x¨+Ω2x=Ω2x−xx˙2−x=G(x,x˙) Say \ddot x + {\Omega ^2}x = {\Omega ^2}x - x{\dot x^2} - x = G(x,\dot x)\;\;\;{\rm{Say}} According to Eq. (4), the direct iterative scheme of Eq. (30) is (32) x¨k+1+Ωk2xk+1=Ωk2xk−xkx˙k2 {\ddot x_{k + 1}} + \Omega _k^2{x_{k + 1}} = \Omega _k^2{x_k} - {x_k}\dot x_k^2 The first approximation x₁(t) and corresponding frequency Ω₀ are obtained from Eq. (32) by substituting k = 0: (33) x1=(A−132(−4+A2)a1,3)cosθ+132(−4+A2)a1,3cos3θ {x_1} = \left({A - {1 \over {32}}(- 4 + {A^2}){a_{1,3}}} \right)\cos \theta + {1 \over {32}}\left({- 4 + {A^2}} \right){a_{1,3}}\cos 3\theta where (34) A1,3=14A3Ω02 {A_{1,3}} = {1 \over 4}{A^3}\Omega _0^2 and (35) Ω0=24−A2. {\Omega _0} = {2 \over {\sqrt {4 - {A^2}}}}. Now, the extended iteration scheme, according to Eq. (3), is as follows: (36) x¨k+1+Ωk2xk+1=G(xk−1,x˙k−1)+Gx(xk−1,x˙k−1)(xk−xk−1)+Gx(xk−1,x˙k−1)(x˙k−x˙k−1); k=1,2,…. {\ddot x_{k + 1}} + \Omega _k^2{x_{k + 1}} = G({x_{k - 1}},{\dot x_{k - 1}}) + {G_x}({x_{k - 1}},{\dot x_{k - 1}})({x_k} - {x_{k - 1}}) + {G_x}({x_{k - 1}},{\dot x_{k - 1}})({\dot x_k} - {\dot x_{k - 1}});\;k = 1,2, \ldots. The second approximation x₂(t) and the corresponding frequency Ω₁ are obtained from Eq. (36) by substituting k = 1: (37) y2=(A+1Ω12(a238+a2524))cosθ−1Ω12(a238cos3θ+a2524cos5θ) {y_2} = \left({A + {1 \over {\Omega _1^2}}\left({{{{a_{23}}} \over 8} + {{{a_{25}}} \over {24}}} \right)} \right)\cos \theta - {1 \over {\Omega _1^2}}\left({{{{a_{23}}} \over 8}\cos 3\theta + {{{a_{25}}} \over {24}}\cos 5\theta} \right) where (38) a23=7A32(4−A62)(−4+A2)−9A58(4−A2)(−4+A2)+A716(4−A2)(−4+A2) −4A3Ω14−A2(−4+A2)+7A5Ω184−A2(−4+A2)+A7Ω1324−A2(−4+A2) A3Ω122(4−A2)(−4+A2)−A5Ω124(4−A2)(−4+A2)+A7Ω1232(4−A2)(−4+A2) \eqalign{& {a_{23}} = {{7{A^3}} \over {2(4 - A62)(- 4 + {A^2})}} - {{9{A^5}} \over {8(4 - {A^2})(- 4 + {A^2})}} + {{{A^7}} \over {16(4 - {A^2})(- 4 + {A^2})}} \cr & \,\,\,\,\,\,\,\,\,\,\,\,\, - {{4{A^3}{\Omega _1}} \over {\sqrt {4 - {A^2}} (- 4 + {A^2})}} + {{7{A^5}{\Omega _1}} \over {8\sqrt {4 - {A^2}} (- 4 + {A^2})}} + {{{A^7}{\Omega _1}} \over {32\sqrt {4 - {A^2}} (- 4 + {A^2})}} \cr & \,\,\,\,\,\,\,\,\,\,\,\,\,{{{A^3}\Omega _1^2} \over {2(4 - {A^2})(- 4 + {A^2})}} - {{{A^5}\Omega _1^2} \over {4(4 - {A^2})(- 4 + {A^2})}} + {{{A^7}\Omega _1^2} \over {32(4 - {A^2})(- 4 + {A^2})}} \cr} (39) a25=A58(4−A2)(−4+A2)−A732(4−A2)(−4+A2)+3A5Ω184−A2(−4+A2)−3A7Ω1324−A2(−4+A2) {a_{25}} = {{{A^5}} \over {8(4 - {A^2})(- 4 + {A^2})}} - {{{A^7}} \over {32(4 - {A^2})(- 4 + {A^2})}} + {{3{A^5}{\Omega _1}} \over {8\sqrt {4 - {A^2}} (- 4 + {A^2})}} - {{3{A^7}{\Omega _1}} \over {32\sqrt {4 - {A^2}} (- 4 + {A^2})}} (40) Ω1=(−128A2+40A4−2A6−2(4096A4−2560A6+528A8−40A10+A12+16384(4−A2) −11264A2(4−A2)+1680A4(4−A2)+32A6(4−A2)−A8(4−A2))1/2)/(2(−1284−A2 +2A24−A2+A44−A2)) \eqalign{& {\Omega _1} = (- 128{A^2} + 40{A^4} - 2{A^6} - 2(4096{A^4} - 2560{A^6} + 528{A^8} - 40{A^{10}} + {A^{12}} + 16384(4 - {A^2}) \cr & \,\,\,\,\,\,\,\,\,\,\,\,\, - 11264{A^2}(4 - {A^2}) + 1680{A^4}(4 - {A^2}) + 32{A^6}(4 - {A^2}) - {A^8}(4 - {A^2}){)^{1/2}})/(2(- 128\sqrt {4 - {A^2}} \cr & \,\,\,\,\,\,\,\,\,\,\,\,\, + 2{A^2}\sqrt {4 - {A^2}} + {A^4}\sqrt {4 - {A^2}})) \cr} In a similar way, the method can be used for higher-order approximations.

An iterative approach to obtain the approximate solution of the ‘quadratic non-linear oscillator’ is presented. The presented technique is very simple for solving algebraic equations analytically, and the approach is different from other existing approaches for adopting truncated Fourier series. This process significantly improves the results.

Here, the first, second and third approximate frequencies Ω₀,Ω₁ and Ω₂ have been calculated, and the results are given in Table 1.

To compare the approximate frequencies, we have also given the existing results determined by Mickens and Ramadhani [9], Belendez et al. [14], Hosen [13] and Haque and Hossain [24], which are shown in Table 2. Fortunately, this current method gives significantly better results than the other formulae.

To show the accuracy, the percentage of errors is calculated by the following definition: Error=|Ωe−ΩkΩe|×100% {\rm{Error}} = \left| {{{{\Omega _e} - {\Omega _k}} \over {{\Omega _e}}}} \right| \times 100\% where Ω_k(k = 0,1,2,...) represents the approximate frequencies obtained by the present method, and Ω_e represents the corresponding exact frequency of the oscillator.

Though the method has been illustrated for a quadratic non-linear oscillator, it is also valid for other oscillators. The same technique has been applied to solve the non-linear jerk oscillator. Herein, we have calculated the approximate frequencies Ω₀, Ω₁ and Ω₂, respectively, in Section 4. The results are given in Table 3 for five particular values of A; to compare the approximate frequencies, we have also given the existing results determined by Gottlieb [?]. We see that the presented technique gives us emphatically better results than the Gottlieb [30] technique.

We know that the basic idea of iterative methods is to construct a sequence of solutions x_k (as well as frequencies Ω_k) that have the property of convergence: xe=limk→∞xk; or ,Ωe=limk→∞Ωk. {x_e} = \mathop {\lim}\limits_{k \to \infty} {x_k};\;{\rm{or}}\;,{\Omega _e} = \mathop {\lim}\limits_{k \to \infty} {\Omega _k}. Here, x_e is the exact solution of the given non-linear oscillator.

In the present method, the solution yields less error in each iterative step compared to the previous iterative step, and finally, |Ω₂ − Ω_e| = |0.91467 − 0.914681 | < ɛ, where ɛ is a small positive number and A is chosen to be unity. From this, it is clear that the adopted method is convergent.

An iterative method of the form represented by Eq. (4), with initial estimate given in Eq. (5) is said to be consistent if limk→∞|xk−xe|=0, or ,limk→∞|Ωk−Ωe|=0 \mathop {\lim}\limits_{k \to \infty} |{x_k} - {x_e}| = 0,\;{\rm{or}}\;,\mathop {\lim}\limits_{k \to \infty} |{\Omega _k} - {\Omega _e}| = 0 In the present analysis, we see that limk→∞|Ωk−Ωe|=0, as |Ω2−Ωe|=0. \mathop {\lim}\limits_{k \to \infty} |{\Omega _k} - {\Omega _e}| = 0,\;{\rm{as}}\;|{\Omega _2} - {\Omega _e}| = 0. Thus, the consistency of the method is achieved.

It is noted that Mickens and Ramadhani [9] found only the second approximate frequencies by the HB method. Belendez et al. [14] found up to the third approximate frequencies by using a modified He's homotopy perturbation method. Again, Hosen [13] found up to the third approximate frequencies by using a modified HB method; Haque and Hossain [24] found up to the fourth approximate frequencies by the iteration method.

In our study, it is seen that the third-order approximate frequency obtained by the adopted method is almost same as the exact frequency. It is found that, in most of the cases, our solution gives significantly better results than other existing results. The advantages of this method include its simplicity and computational efficiency.