An Optimal Iteration Method for Strongly Nonlinear Oscillators

Abstract

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 [1–3]. 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 [5–8], 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 𝑢₀(τ) can be chosen in the general form 𝑢₀(𝜏)=^𝑚_𝑖=1𝐶_𝑖𝑓_𝑖(𝜏),(2.4) where C_i 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, ((iii))𝑢𝑢^𝑛+1_𝑛+1(𝜏)+𝐹_𝑛𝜏,𝐶₁,𝐶₂,…,𝐶_𝑚+𝐶′=0,(𝜏)+𝐺_𝑛𝜏,𝐶₁,𝐶₂,…,𝐶_𝑚+𝐶^𝜏+𝐶^=0,(2.5) where 𝐺𝐹_𝑛𝜏,𝐶₁,𝐶₂,…,𝐶_𝑚=𝑓𝑢_𝑛(𝜏),𝑢^𝑛(𝜏),𝑢^𝑛(𝜏)𝑑𝜏,_𝑛𝜏,𝐶₁,𝐶₂,…,𝐶_𝑚=𝐹_𝑛𝜏,𝐶₁,𝐶₂,…,𝐶_𝑚𝑑𝜏.(2.6)

From the initial conditions (2.2), we consider ((iii))𝐹𝐺_𝑛𝑛0,𝐶₁,𝐶₂,…,𝐶_𝑚=0,0,𝐶₁,𝐶₂,…,𝐶_𝑚=−𝐴(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(𝜏)=−𝐺_𝑛𝜏,𝐶₁,𝐶₂,…,𝐶_𝑚,(2.8) where the constants 𝐶₁,𝐶₂,…,𝐶_𝑚 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 𝐽𝐶₁,𝐶₂,…,𝐶_𝑚=^𝑇0𝑢^𝑛+𝑓𝑢_𝑛,𝑢^𝑛,𝑢^𝑛²𝑑𝜏(2.9) is minimum, one can obtain the optimal values of the unknown constants. Taking into consideration (2.7), the constants C_i, 𝑖=1,2,…,𝑚 can be determined in this case from the equations (conditioned minimum) 𝜕𝐽𝜕𝐶_𝑗+𝜆₁𝜕𝐹_𝑛0,𝐶₁,𝐶₂,…,𝐶_𝑚𝜕𝐶_𝑗+𝜆₂𝜕𝐺_𝑛0,𝐶₁,𝐶₂,…,𝐶_𝑚𝜕𝐶_𝑗=0,𝑗=3,4,…,𝑚,(2.10) where 𝜆₁=𝜕𝐽/𝜕𝐶¹𝜕𝐺_𝑛/𝜕𝐶₂−𝜕𝐽/𝜕𝐶₂𝜕𝐺_𝑛/𝜕𝐶₁𝜆₂𝜕𝐺_𝑛/𝜕𝐶₁𝜕𝐹_𝑛/𝜕𝐶₂−𝜕𝐺_𝑛/𝜕𝐶₂𝜕𝐹_𝑛/𝜕𝐶₁,=𝜕𝐽/𝜕𝐶¹𝜕𝐹_𝑛/𝜕𝐶₂−𝜕𝐽/𝜕𝐶₂𝜕𝐹_𝑛/𝜕𝐶₁𝜕𝐹_𝑛/𝜕𝐶₁𝜕𝐺_𝑛/𝜕𝐶₂−𝜕𝐹_𝑛/𝜕𝐶₂𝜕𝐺_𝑛/𝜕𝐶₁(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: 𝜕𝐽𝜕𝐶_𝑗−𝜕𝐽/𝜕𝐶¹𝜕𝐺_𝑛/𝜕𝐶₁𝜕𝐺_𝑛𝜕𝐶_𝑗=0,𝑗=2,3,…,𝑚.(2.12)

Therefore, the solution (2.8) with the known constants 𝐶₁,𝐶₂,…,𝐶_𝑚 is well determined.

In the present paper we consider a nonlinear oscillator with cubic and harmonic restoring force 𝑢̈+𝑢+𝑎𝑢³+𝑏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Ω²𝑢+𝑎𝑢³+𝑏sin𝑢.(2.17)

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

We consider the initial approximation in the form 𝑢₀(𝜏)=𝐶₁cos𝜏+2𝐶₂cos3𝜏+2𝐶₃cos5𝜏+2𝐶₄cos7𝜏,(2.19) where C₁, C₂, C₃, and C₄ are unknown constants at this moment.

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

The first term in the right-hand side of (2.23) becomes 𝑓𝐶₁cos𝜏=−1Ω²𝐶¹cos𝜏+𝑎𝐶³¹4(cos3𝜏+3cos𝜏)+𝑏sin𝐶¹cos𝜏.(2.24)

The last term in (2.24) can be expanded in the power series sin𝐶₁cos𝜏=𝐶₁cos𝜏−13!𝐶³¹cos³𝜏+15!𝐶⁵¹cos⁵𝜏−17!𝐶⁷¹cos⁷𝜏+19!𝐶⁹¹cos⁹𝜏+⋯.(2.25)

Substituting (2.25) into (2.24), after some simple manipulations we obtain 𝑓𝐶₁cos𝜏=𝛼₁cos𝜏+𝛼₃cos3𝜏+𝛼₅cos5𝜏+𝛼₇cos7𝜏+𝛼₉cos9𝜏+⋯,(2.26) where 𝛼₁=−𝐶¹Ω²1+34𝑎𝐶²¹+𝑏1−𝐶²¹𝛼8+1𝐶9⁴¹2−9𝐶21⁶¹6+73𝐶72⁸¹80+⋯₃;=−𝐶³¹𝛼Ω²14𝑎−𝑏214−3𝐶8¹4+15𝐶3⁴¹60−110𝐶5⁶¹920+⋯₅;=−𝑏𝐶⁵¹𝛼Ω²19120−46𝐶0²¹80+258𝐶0⁴¹480+⋯₇;=𝑏𝐶⁷¹Ω²3221560−103𝐶21²¹920+⋯;𝛼⁹=−𝑏𝐶⁹¹Ω²928917280+⋯.(2.27)

The last term in the right-side of (2.23) is 𝑓_𝑢𝐶₁cos𝜏=−1Ω²1+3𝑎𝐶²¹cos²𝜏+𝑏cos𝐶₁cos𝜏.(2.28)

In (2.28), the last term can be written as cos𝐶₁cos𝜏=1−𝐶²¹cos²𝜏2!+𝐶⁴¹cos⁴𝜏4!−𝐶⁶¹cos⁶𝜏6!+𝐶⁸¹cos⁸𝜏8!+⋯.(2.29)

Substituting (2.29) into (2.28), we obtain 𝑓_𝑢𝐶₁cos𝜏=𝛽₀+𝛽₂cos2𝜏+𝛽₄cos4𝜏+𝛽₆cos6𝜏+𝛽₈cos8𝜏+⋯,(2.30) where 𝛽₀=−1Ω²1+32𝑎𝐶²¹+𝑏1−𝐶²¹𝛽4+𝐶64⁴¹−2𝐶30⁶¹4+14𝐶74⁸¹56+⋯₂;=1Ω²32𝑎𝐶²¹−1𝛽4𝐶²¹1−𝐶12²¹+3𝐶8⁴¹4−23𝐶0⁶¹40+⋯₄;=𝑏𝐶⁴¹𝛽192Ω²1−𝐶20²¹+9𝐶6⁴¹0+⋯₆;=−𝑏𝐶⁶¹2304Ω²1−𝐶28²¹+⋯;𝛽⁸=𝑏𝐶⁸¹5160960Ω²(1+⋯).(2.31)

Substituting (2.24) and (2.30) into (2.23), we obtain the expression 𝑓𝑢₀=+𝛼₁𝛼+𝛽₂+𝛽₄𝐶₂+𝛽₄+𝛽₆𝐶₃+𝛽₆+𝛽₈𝐶₄cos𝜏₃+𝛼+2𝛽₀+𝛽₆𝐶₂+𝛽₂+𝛽₈𝐶₃+𝛽₄𝐶₄cos3𝜏₅+𝛼+𝛽₂+𝛽₈𝐶₂+2𝛽₀𝐶₃+𝛽₂𝐶₄cos5𝜏₇+𝛼+2𝛽₄𝐶₂+2𝛽₂𝐶₃+2𝛽₀𝐶₄cos7𝜏₉+2𝛽₆𝐶₂+2𝛽₄𝐶₃+2𝛽₂𝐶₄cos9𝜏+⋯.(2.32) Equation (2.5)(i) becomes 𝑢^1(𝜏)=−−1𝛼₁+𝛽₂+𝛽₄𝐶₂+𝛽₄+𝛽₆𝐶₃+𝛽₆+𝛽₈𝐶₄sin𝜏3𝛼³−1+2𝛽₀+𝛽₆𝐶₂+𝛽₂+𝛽₈𝐶₃+𝛽₄𝐶₄sin3𝜏5𝛼⁵−1+𝛽₂+𝛽₈𝐶₂+2𝛽₀𝐶₃+𝛽₂𝐶₄sin5𝜏7𝛼⁷−1+2𝛽₄𝐶₂+2𝛽₂𝐶₃+2𝛽₀𝐶₄sin7𝜏9𝛼⁹+2𝛽₆𝐶₂+2𝛽₄𝐶₃+2𝛽₂𝐶₄sin9𝜏+⋯.(2.33)

Finally, (2.8) becomes 𝑢₁(𝜏)=+𝛼1₁+𝛽₂+𝛽₄𝐶₂+𝛽₄+𝛽₆𝐶₃+𝛽₆+𝛽₈𝐶₄cos𝜏9𝛼³+1+2𝛽₀+𝛽₆𝐶₂+𝛽₂+𝛽₈𝐶₃+𝛽₄𝐶₄cos3𝜏25𝛼⁵+1+𝛽₂+𝛽₈𝐶₂+2𝛽₀𝐶₃+𝛽₂𝐶₄cos5𝜏49𝛼⁷+1+2𝛽₄𝐶₂+2𝛽₂𝐶₃+2𝛽₀𝐶₄cos7𝜏81𝛼⁹+2𝛽₆𝐶₂+2𝛽₄𝐶₄+2𝛽₂𝐶₄cos9𝜏.(2.34) From (2.33) we obtain that (2.7)(i) becomes identity and (2.7)(ii) becomes 𝛼₁+19𝛼³+125𝛼⁵+149𝛼⁷+181𝛼⁹+𝐶₂29𝛽⁰+2625𝛽²+5149𝛽⁴+1181𝛽⁶+1+𝐶₃25𝛽⁸225𝛽⁰+46471𝛽²+8831𝛽⁴+𝛽₆+19𝛽⁸+𝐶⁴429𝛽⁰+2103215𝛽²+19𝛽⁴+𝛽₆+𝛽₈−𝐴=0.(2.35)

The frequency Ω and the constants 𝐶₁, 𝐶₂, 𝐶₃, and 𝐶₄ 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 1–6 for several cases.

Figure 1 

Comparison between the approximate solution (2.34) and numerical solution of (2.13) in Case a: 𝑎=𝑏=𝐴=1: dashed red line: numerical solution, dashed blue line approximate solution.

Figure 2 

Comparison between the approximate solution (2.34) and numerical results of (2.13) in terms of phase plane in Case a: 𝑎=𝑏=𝐴=1: dashed red line: numerical solution, dashed blue line approximate solution.

Figure 3 

The error between the numerical and approximate solution (2.34) in Case a: 𝑎=𝑏=𝐴=1.

Figure 4 

Comparison between the approximate solution (2.34) and numerical solution of (2.13) in Case b: 𝑎=𝑏=1,𝐴=2: dashed red line: numerical solution, dashed blue line approximate solution.

Figure 5 

Comparison between the approximate solution (2.34) and numerical results of (2.13) in terms of phase plane in Case b: 𝑎=𝑏=1,𝐴=2: dashed red line: numerical solution, dashed blue line approximate solution.

Figure 6 

The error between the numerical and approximate analytical solution (2.34) in case b: 𝑎=𝑏=1,𝐴=2.

Case a. For 𝑎=1, 𝑏=1, 𝐴=1, following the procedure described above we obtain the approximate periodic solution of (2.13) in the form 𝑢₁(𝑡)=+0.09.8080309040539979c4ocsoΩs7𝑡Ω+𝑡0−.001.10301000234518c1ocso3sΩ9Ω𝑡+𝑡,0.000326978cos5Ω𝑡(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: 𝐸𝑟(𝑡)=𝑢_𝑁(𝑡)−𝑢₁(𝑡),(3.2) where 𝑢_𝑁(t) is the numerical result and 𝑢₁(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.09.4070005224311321c2ocsoΩs7𝑡Ω+𝑡0−.005.20101071922736c3o5sc3oΩs9𝑡Ω+𝑡0,.001198712cos5Ω𝑡(3.3) where Ω=2.12453. Comparisons between the approximate and numerical results for Case b are presented in Figures 4–6.It can be seen from Figures 1–6 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.