We introduce a new method, namely, the Optimal Iteration Perturbation Method (OIPM), to solve nonlinear differential equations of oscillators with cubic and harmonic restoring force. We illustrate that OIPM is very effective and convenient and does not require linearization or small perturbation. Contrary to conventional methods, in OIPM, only one iteration leads to high accuracy of the solutions. The main advantage of this approach consists in that it provides a convenient way to control the convergence of approximate solutions in a very rigorous way and allows adjustment of convergence regions where necessary. A very good agreement was found between approximate and numerical solutions, which prove that OIPM is very efficient and accurate.

1. Introduction

Mathematical modelling of many physical systems leads to nonlinear ordinary or partial differential equations in various fields of physics, mathematics, or engineering. An effective method is required to analyze the mathematical model which provides solutions conforming to physical reality. In many cases, it is possible to replace a nonlinear differential equation by a corresponding linear differential equation that approximates closely the original one to give useful results. In general, the study of nonlinear differential equations is restricted to a variety of special classes of equations and the method of solution usually involves one or more techniques to achieve analytical approximations to the solutions. Solving the governing equations of nonlinear oscillators has been one of the most time-consuming and difficult affairs among researchers. Therefore, many researchers and scientists of both vibrations and mathematics have recently paid much attention to find and develop approximate solutions. Perturbation methods are well established tools to study diverse aspects of nonlinear problems [13]. However, the use of perturbation theory in many important practical problems is invalid, or it simply breaks down for parameters beyond a certain specified range. Therefore, new analytical techniques should be developed to overcome these shortcomings. Such a new technique should work over a larger range of parameters and yield accurate analytical approximate solutions beyond the coverage and ability of the classical perturbation methods.

It is noted that several methods have been used to obtain approximate solutions for strongly nonlinear oscillators. An interesting approach which combines the harmonic balance method and linearization of nonlinear oscillation equation was proposed in [4]. There also exists a wide range of literature dealing with approximate periodic solutions for nonlinear problems with large parameters by using a mixture of methodologies: the variational iteration method [58], some linearization methods [9, 10], the optimal homotopy asymptotic method [11], the optimal parametric iteration method [12], some modified Lindstedt-Poincare methods [13, 14], or a simple approach [15].

In this paper, coupling the iteration perturbation method [16] with the least square technology, a new approach, namely, the Optimal Iteration Perturbation Method (OIPM), is proposed to find explicit analytical periodic solutions to nonlinear oscillators with cubic and harmonic restoring force. Recently, in the same way, the variational iteration method [5] and the homotopy perturbation method [17] have been coupled with the least square technology resulting in two new powerful methods, namely, the optimal variational iteration method (OVIM) [7] and the optimal homotopy perturbation method (OHPM) [18].

The efficiency of the present procedure is proved while an accurate solution is explicitly analytically obtained in an iterative way after only one iteration. The proposed method does not require a small parameter into the equation and provides a convenient and rigorous way to optimally control the convergence of the solutions by means of a finite number of unknown parameters.

2. Formulation and Solution Approach

In this work, we consider a nonlinear oscillator in the form 𝑢+𝑓𝑢,𝑢,𝑢=0,(2.1) with initial conditions 𝑢(0)=𝐴,𝑢(0)=0,(2.2) where prime denotes derivative with respect to variable 𝜏.

For (2.1) and (2.2) we propose the following iteration scheme: 𝑢𝑛+1𝑢+𝑓𝑛,𝑢𝑛,𝑢𝑛=0,𝑛=0,1,2,,(2.3) where the initial approximation 𝑢0(τ) can be chosen in the general form 𝑢0(𝜏)=𝑚𝑖=1𝐶𝑖𝑓𝑖(𝜏),(2.4) where Ci are unknown constants, m is a positive integer number, and the functions 𝑓𝑖 are trigonometric functions sine or/and cosine in case of nonlinear oscillators.

Integrating (2.3) twice with respect to τ, we have, respectively, (i)𝑢𝑛+1(𝜏)+𝐹𝑛𝜏,𝐶1,𝐶2,,𝐶𝑚+𝐶=0,(ii)𝑢𝑛+1(𝜏)+𝐺𝑛𝜏,𝐶1,𝐶2,,𝐶𝑚+𝐶𝜏+𝐶=0,(2.5) where 𝐹𝑛𝜏,𝐶1,𝐶2,,𝐶𝑚=𝑓𝑢𝑛(𝜏),𝑢𝑛(𝜏),𝑢𝑛(𝐺𝜏)𝑑𝜏,𝑛𝜏,𝐶1,𝐶2,,𝐶𝑚=𝐹𝑛𝜏,𝐶1,𝐶2,,𝐶𝑚𝑑𝜏.(2.6)

From the initial conditions (2.2), we consider (i)𝐹𝑛0,𝐶1,𝐶2,,𝐶𝑚(=0,ii)𝐺𝑛0,𝐶1,𝐶2,,𝐶𝑚=𝐴(2.7) such that the integration constants 𝐶 and 𝐶 into (2.7)(i) and (2.5)(ii) become 𝐶=𝐶=0.

In this way, the approximate solution of n+1 order can be written in the form 𝑢𝑛+1(𝜏)=𝐺𝑛𝜏,𝐶1,𝐶2,,𝐶𝑚,(2.8) where the constants 𝐶1,𝐶2,,𝐶𝑚 which are considered in the initial approximation (2.4) can be identified via various methods, such as, for example, the least square method, the Galerkin method, the Ritz method, and the collocation method. For example, imposing that the residual functional given by 𝐽𝐶1,𝐶2,,𝐶𝑚=𝑇0𝑢𝑛𝑢+𝑓𝑛,𝑢𝑛,𝑢𝑛2𝑑𝜏(2.9) is minimum, one can obtain the optimal values of the unknown constants. Taking into consideration (2.7), the constants Ci, 𝑖=1,2,,𝑚 can be determined in this case from the equations (conditioned minimum) 𝜕𝐽𝜕𝐶𝑗+𝜆1𝜕𝐹𝑛0,𝐶1,𝐶2,,𝐶𝑚𝜕𝐶𝑗+𝜆2𝜕𝐺𝑛0,𝐶1,𝐶2,,𝐶𝑚𝜕𝐶𝑗=0,𝑗=3,4,,𝑚,(2.10) where 𝜆1=𝜕𝐽/𝜕𝐶1𝜕𝐺𝑛/𝜕𝐶2𝜕𝐽/𝜕𝐶2𝜕𝐺𝑛/𝜕𝐶1𝜕𝐺𝑛/𝜕𝐶1𝜕𝐹𝑛/𝜕𝐶2𝜕𝐺𝑛/𝜕𝐶2𝜕𝐹𝑛/𝜕𝐶1,𝜆2=𝜕𝐽/𝜕𝐶1𝜕𝐹𝑛/𝜕𝐶2𝜕𝐽/𝜕𝐶2𝜕𝐹𝑛/𝜕𝐶1𝜕𝐹𝑛/𝜕𝐶1𝜕𝐺𝑛/𝜕𝐶2𝜕𝐹𝑛/𝜕𝐶2𝜕𝐺𝑛/𝜕𝐶1(2.11) and if (2.7)(i) is not identity. Now, if (2.7)(i) becomes identity, the constants 𝐶𝑖,𝑖=1,2,,𝑚 then can be determined from (2.7)(ii) and from the following equations: 𝜕𝐽𝜕𝐶𝑗𝜕𝐽/𝜕𝐶1𝜕𝐺𝑛/𝜕𝐶1𝜕𝐺𝑛𝜕𝐶𝑗=0,𝑗=2,3,,𝑚.(2.12)

Therefore, the solution (2.8) with the known constants 𝐶1,𝐶2,,𝐶𝑚 is well determined.

In the present paper we consider a nonlinear oscillator with cubic and harmonic restoring force ̈𝑢+𝑢+𝑎𝑢3+𝑏sin𝑢=0,(2.13) where a and b are known constants and dot denotes derivative with respect to time t. The initial conditions are given by 𝑢(0)=𝐴,̇𝑢(0)=0.(2.14)

If Ω is the frequency of the system described by (2.13) and introducing a new independent variable 𝜏=Ω𝑡(2.15) then (2.13) becomes 𝑢+𝑓(𝑢)=0,(2.16) where =𝑑/𝑑𝜏 and 1𝑓(𝑢)=Ω2𝑢+𝑎𝑢3+𝑏sin𝑢.(2.17)

The initial conditions (2.14) become 𝑢(0)=𝐴,𝑢(0)=0.(2.18)

We consider the initial approximation in the form 𝑢0(𝜏)=𝐶1cos𝜏+2𝐶2cos3𝜏+2𝐶3cos5𝜏+2𝐶4cos7𝜏,(2.19) where C1, C2, C3, and C4 are unknown constants at this moment.

For 𝑛=0 into (2.3) we obtain the first iteration given by 𝑢1𝑢+𝑓0=0(2.20) but it is difficult to calculate 𝑓(𝑢0) with u0 given by (2.19). Now, the function f can be expanded in a series using the well-known formula 𝑓𝑡0𝑡+=𝑓0+𝑓1!𝑢𝑡0+,(2.21) where 𝑓𝑢=𝑑𝑓/𝑑𝑢. In the following, we consider 𝑡0=𝐶1cos𝜏,=2𝐶2cos3𝜏+2𝐶3cos5𝜏+2𝐶4cos7𝜏(2.22) such that, from (2.19), (2.21), and (2.22), we obtain 𝑓𝑢0𝐶=𝑓1+cos𝜏2𝐶2cos3𝜏+2𝐶3cos5𝜏+2𝐶4𝑓cos7𝜏𝑢𝐶1cos𝜏.(2.23)

The first term in the right-hand side of (2.23) becomes 𝑓𝐶11cos𝜏=Ω2𝐶1cos𝜏+𝑎𝐶314𝐶(cos3𝜏+3cos𝜏)+𝑏sin1cos𝜏.(2.24)

The last term in (2.24) can be expanded in the power series 𝐶sin1cos𝜏=𝐶11cos𝜏𝐶3!31cos31𝜏+𝐶5!51cos51𝜏𝐶7!71cos71𝜏+𝐶9!91cos9𝜏+.(2.25)

Substituting (2.25) into (2.24), after some simple manipulations we obtain 𝑓𝐶1cos𝜏=𝛼1cos𝜏+𝛼3cos3𝜏+𝛼5cos5𝜏+𝛼7cos7𝜏+𝛼9cos9𝜏+,(2.26) where 𝛼1𝐶=1Ω231+4𝑎𝐶21𝐶+𝑏1218+𝐶41𝐶19261+𝐶921681;𝛼737280+3𝐶=31Ω2141𝑎𝑏𝐶241+𝐶38441𝐶1536061;𝛼1105920+5=𝑏𝐶51Ω21𝐶192021+𝐶4608041;𝛼2580480+7=𝑏𝐶71Ω21𝐶3225602110321920+;𝛼9=𝑏𝐶91Ω21.92897280+(2.27)

The last term in the right-side of (2.23) is 𝑓𝑢𝐶11cos𝜏=Ω21+3𝑎𝐶21cos2𝐶𝜏+𝑏cos1cos𝜏.(2.28)

In (2.28), the last term can be written as 𝐶cos1𝐶cos𝜏=121cos2𝜏+𝐶2!41cos4𝜏𝐶4!61cos6𝜏+𝐶6!81cos8𝜏8!+.(2.29)

Substituting (2.29) into (2.28), we obtain 𝑓𝑢𝐶1cos𝜏=𝛽0+𝛽2cos2𝜏+𝛽4cos4𝜏+𝛽6cos6𝜏+𝛽8cos8𝜏+,(2.30) where 𝛽01=Ω231+2𝑎𝐶21𝐶+𝑏1214+𝐶41𝐶6461+𝐶230481;𝛽147456+2=1Ω232𝑎𝐶2114𝐶21𝐶121+𝐶1241𝐶38461;𝛽23040+4=𝑏𝐶41192Ω2𝐶121+𝐶2041;𝛽960+6=𝑏𝐶612304Ω2𝐶12128+;𝛽8=𝑏𝐶815160960Ω2(1+).(2.31)

Substituting (2.24) and (2.30) into (2.23), we obtain the expression 𝑓𝑢0=𝛼1+𝛽2+𝛽4𝐶2+𝛽4+𝛽6𝐶3+𝛽6+𝛽8𝐶4+𝛼cos𝜏3+2𝛽0+𝛽6𝐶2+𝛽2+𝛽8𝐶3+𝛽4𝐶4+𝛼cos3𝜏5+𝛽2+𝛽8𝐶2+2𝛽0𝐶3+𝛽2𝐶4+𝛼cos5𝜏7+2𝛽4𝐶2+2𝛽2𝐶3+2𝛽0𝐶4+𝛼cos7𝜏9+2𝛽6𝐶2+2𝛽4𝐶3+2𝛽2𝐶4cos9𝜏+.(2.32) Equation (2.5)(i) becomes 𝑢1𝛼(𝜏)=1+𝛽2+𝛽4𝐶2+𝛽4+𝛽6𝐶3+𝛽6+𝛽8𝐶41sin𝜏3𝛼3+2𝛽0+𝛽6𝐶2+𝛽2+𝛽8𝐶3+𝛽4𝐶41sin3𝜏5𝛼5+𝛽2+𝛽8𝐶2+2𝛽0𝐶3+𝛽2𝐶41sin5𝜏7𝛼7+2𝛽4𝐶2+2𝛽2𝐶3+2𝛽0𝐶41sin7𝜏9𝛼9+2𝛽6𝐶2+2𝛽4𝐶3+2𝛽2𝐶4sin9𝜏+.(2.33)

Finally, (2.8) becomes 𝑢1𝛼(𝜏)=1+𝛽2+𝛽4𝐶2+𝛽4+𝛽6𝐶3+𝛽6+𝛽8𝐶4+1cos𝜏9𝛼3+2𝛽0+𝛽6𝐶2+𝛽2+𝛽8𝐶3+𝛽4𝐶4+1cos3𝜏𝛼255+𝛽2+𝛽8𝐶2+2𝛽0𝐶3+𝛽2𝐶4+1cos5𝜏𝛼497+2𝛽4𝐶2+2𝛽2𝐶3+2𝛽0𝐶4+1cos7𝜏𝛼819+2𝛽6𝐶2+2𝛽4𝐶4+2𝛽2𝐶4cos9𝜏.(2.34) From (2.33) we obtain that (2.7)(i) becomes identity and (2.7)(ii) becomes 𝛼1+19𝛼3+1𝛼255+1𝛼497+1𝛼819+𝐶229𝛽0+26𝛽252+51𝛽494+11𝛽816+1𝛽258+𝐶32𝛽250+67𝛽4412+83𝛽814+𝛽6+19𝛽8+𝐶42𝛽490+131𝛽20252+19𝛽4+𝛽6+𝛽8𝐴=0.(2.35)

The frequency Ω and the constants 𝐶1, 𝐶2, 𝐶3, and 𝐶4 are determined by means of a collocation-type method.

3. Numerical Examples

We will illustrate the applicability, accuracy, and effectiveness of the proposed approach by comparing the analytical approximate periodic solution with numerical integration results obtained using a fourth-order Runge-Kutta method. The comparison is made in terms of displacements and phase plane. The error of the solution has been also computed. The results of these comparisons are presented in Figures 16 for several cases.

Case a. For 𝑎=1, 𝑏=1, 𝐴=1, following the procedure described above we obtain the approximate periodic solution of (2.13) in the form 𝑢1(𝑡)=0.988394597cosΩ𝑡+0.011310241cos3Ω𝑡+0.000326978cos5Ω𝑡+0.000003994cos7Ω𝑡0.00003581cos9Ω𝑡,(3.1) where Ω=1.61923. In Figure 1 is presented a comparison between the approximate solution (3.1) and the solution obtained through numerical simulations. Moreover, Figure 2 presents a comparison between the approximate solution (3.1) and the numerical results in terms of phase plane. In order to provide a comprehensive evidence of the accuracy of the results, the error of the solution has been computed: 𝐸𝑟(𝑡)=𝑢𝑁(𝑡)𝑢1(𝑡),(3.2) where 𝑢𝑁(t) is the numerical result and 𝑢1(t) is the approximate solution given by (2.34). A graphical representation of the error in the Case a is presented in Figure 3.

Case b. For 𝑎=1,𝑏=1,𝐴=2, following the same procedure we obtain 𝑢1(𝑡)=1.947052312cosΩ𝑡+0.052117923cos3Ω𝑡+0.001198712cos5Ω𝑡0.000241312cos7Ω𝑡0.000127635cos9Ω𝑡,(3.3) where Ω=2.12453. Comparisons between the approximate and numerical results for Case b are presented in Figures 46.It can be seen from Figures 16 that the results obtained using OIPM are almost identical with those obtained through numerical simulations.

4. Conclusions

In this paper we have developed an analytical treatment of strongly nonlinear oscillators with cubic and harmonic restoring force using a new approximate analytical technique, namely, the Optimal Iteration Perturbation Method (OIPM). This method accelerates the convergence of the solutions since after only one iteration we achieved very accurate results. The proposed approach is an iterative procedure, and iterations are preformed in a very simple manner by identifying optimally some coefficients and therefore very good approximations are obtained in few terms. Actually, the capital strength of OIPM is its fast convergence. An excellent agreement of the approximate periodic solutions and frequencies with the exact ones has been demonstrated. Two examples are given, and the results reveal that our procedure is very effective, simple, and accurate. This paper demonstrates the general validity and the great potential of the OIPM for solving strongly nonlinear problems.