Nonlinear oscillators dynamics using optimal and modified homotopy perturbation method

Abstract

This study investigates the several strongly nonlinear oscillators (Van der Pol, Duffing and Rayleigh) for which we have applied the most powerful and advanced optimal semi-analytical technique: optimal and modified homotopy perturbation method (OM-HPM) for the convergent semi-analytical series solution. The numerical simulation demonstrates the high accuracy of the OM-HPM, which is straightforward, does not require any domain decomposition, special transformation, or pade approximations to get the convergent series solution. The key features for the high accuracy of the OM-HPM lies on the best optimal auxiliary linear operator. Therefore, OM-HPM offering a valuable tool for engineers and researcher to analyze the complex nonlinear oscillator.

Introduction

Many scientific fields, including physics, engineering, and structures including thin-walled structures ¹ and active structure ² use nonlinear oscillators. Many nonlinear differential equations appear in these fields and studying these equations is a major issue. Lots of efforts have been made to explore them in literature, and numerous analytical and numerical techniques have been proposed and employed to solve such problems.^3–7 In literature, several analytical methods have been utilized since an exact solution could be too complex to be applied in real-world situations. For instance, the Hamiltonian approach,^8,9 the variational iteration method,^10–12 the homotopy perturbation method,^13–17 the variational approach,^18–21 in addition to many other methods.^22–25

Nonlinear oscillators come in two varieties: weakly and extremely nonlinear. Weak nonlinearity is indicated by a small coefficient of the nonlinear term; however, some authors have claimed that this is not always enough and that the amplitude and starting conditions must also be taken into consideration. Weak nonlinearity can be quantified by treating it as a disturbance of the corresponding linear oscillator. There are several perturbation techniques available, but not all of them yield consistently reliable results. Modifications have been developed to alleviate these shortcomings.^26,27

Numerous researchers have shown interest in a number of notable oscillators that have important applications in physics and engineering.^28–32 For example, the excited spring pendulum has been studied in detail utilizing the so-called multi-scale approach, which enables the authors to estimate the equations of motion for the system. ²⁹ Additionally, the Van der Poland Duffing oscillators have been considered in a number of papers.^29–31 It is clear from these references that several perturbation approaches were investigated in order to approximate the solution to the dynamic equations of the Van der Pol oscillator.^29,30 On the other hand, the Duffing oscillator was also studied in References 29–32, where the authors used the effective analytical approach known as the fourth-order Runge–Kutta method to determine the solution and compare it with a numerical solution.

Semi-analytical approaches, which combine analytical and numerical techniques, are becoming more and more common for approximating and solving nonlinear problems. They are useful for comprehending and handling complicated nonlinear occurrences in a variety of domains because they provide a balance between accuracy and computing efficiency. Among such semi-analytical approaches, homotopy perturbation method (HPM)^13,14,32 is one of the most powerful for solving the highly nonlinear oscillators and is widely use to solve such problems. However, the accuracy and convergence of the HPM series solution are slower for few problems.^33,34 Even sometimes HPM^33,34 fails to provide us the convergent series solution for few highly nonlinear problems. Moreover, others semi-analytical methods including variational iteration method (VIM) ³⁵ and differential transformation method (DTM)^34,35 are also fails to provide convergent series solutions. In recent literature, there is most powerful and advance technique: an optimal and modified homotopy perturbation method (OM-HPM). ³⁶ Which is nothing but the optimal form of the classical HPM^13,14,32 and it is well established for highly nonlinear differential equations to get the convergence series solution. However, it is seen that the accuracy and convergence of the HPM series solution are slower for few problems.^33,34 Even sometimes HPM^33,34 fails to provide us the convergent series solution for some of the highly nonlinear problems, particularly, for the van der pol oscillator. Several researchers^11,37–41 put effort to solve such van der pol oscillator using other contemporary semi-analytical methods including variational iteration method (VIM), ³⁵ optimal VIM, ⁴² differential transformation method (DTM).^37,38 Unfortunately, they failed to provide convergent series solutions. In the literature, we have added an advanced semi-analytical method, namely, Optimal and Modified Homotopy Perturbation Method (OM-HPM). ³⁶ The traditional Homotopy Perturbation Method (HPM)^{3,13–16,19,24–29} is updated by defining the linear operator as an auxiliary linear operator, and then optimized this operator by minimizing the residual error. Unlike traditional methods, OM-HPM ³⁶ does not require any auxiliary parameters or functions. In addition, it can be directly applied to singular or non-singular highly nonlinear ordinary differential equations without decomposition, special transformations, or pade approximation. Therefore, in this study, our aim is to analyze the dynamics of some of the strong nonlinear oscillators, namely, duffing oscillator (DO), Rayleigh oscillator (RO), and van der Pol oscillator (vdPO), using this advance optimal analytical technique.

The rest of this paper is organized as follows: In section 2, we present the methodology of the method used. In section 3, analytical and numerical methods have been applied to solve three oscillators (Duffing, Van der Pol and Rayleigh). Finally, we close our paper by a conclusion.

Methodology

In order to solve the considered nonlinear problems, we adopted the recently proposed optimal and modified homotopy perturbation method (OM-HPM). ³⁶ The method is follows:

For the nonlinear differential equation of the form

$\mathcal{N}\left[x\right]=\mathcal{A}\left[x\right]+g\left(t\right)=0,t\in \mathrm{\Omega }$

(1)with the boundary condition

$\mathcal{B}\left(x,\frac{\partial x}{\partial t}\right)=0,t\in \mathrm{\Gamma }$

(2)where $\mathcal{N}$ , nonlinear operator; $\mathcal{A}$ , general differential operator; $\mathcal{B}$ , boundary operator; $\mathrm{\Gamma }$ , boundary of the domain $\mathrm{\Omega }$ ; and $g$ and $x$ are the known analytic and unknown functions, respectively. The zero-th order homotopy equation as

$\left(1-p\right)\hspace{0.17em}\mathcal{L}\left[F\left(t,p\right)-x_0\left(t\right)\right]+p\hspace{0.17em}\mathcal{N}\left[F\left(t,p\right)\right]=0$

(3)where $p\in \left[0,1\right]$ , embedding parameter; $F\left(t,p\right)$ , unknown; $x_0\left(t\right)$ , initial approximation; and $\mathcal{L}$ , auxiliary linear operator of the form

$\begin{array}{l}\mathcal{L}\left[\mathrm{x}\left(t\right)\right]=\frac{d^{n}x}{d{t}^n}+a_{n-1}\frac{d^{n-1}x}{d{t}^{n-1}}+a_{n-2}\frac{d^{n-2}x}{d{t}^{n-2}}+......+a_1\frac{dx}{dt}+{\left(-1\right)}^n{a}_{0}x\left(t\right)\\ \hspace{0.17em}\text{with}\\ \hspace{0.17em}{a}_{n-k}=\frac{{\left(-1\right)}^k}{k!}{B}_k\left(s_1,-1!s_2,2!s_3,.......,{\left(-1\right)}^{k-1}\left(k-1\right)!s_k\right),\end{array}$

(4)where $B_k$ is the Bell’s polynomial, ⁸ $s_k={{\displaystyle \sum }}_{i=1}^n\hspace{0.17em}{{\lambda }_i}^k$ and $a_0={\displaystyle \prod }_{i=1}^n\hspace{0.17em}{\lambda }_i$ , ${\lambda }_i$ ’s are the roots of the auxiliary equation of $\mathcal{L}\left[x\left(t\right)\right]=0$ . The auxiliary linear operator $\mathcal{L}$ satisfied the property $\mathcal{L}\left[\mathrm{x}\left(\mathrm{t}\right)\right]=0$ for $\mathrm{x}\left(\mathrm{t}\right)=0$ . For $p=0$ , we get from equation (3):

$F\left(t,0\right)=x_0\left(t\right).$

Also when $p=1$ from equation (3) we have

$F\left(t,1\right)=x\left(t\right),$

this represents the solution of the given equation (1). As $p$ increase from $0$ to $1$ , the initial solution $x_0\left(t\right)$ continuously deforms to the final solution $x\left(t\right)$ of the problem. That is, the Taylor series

$F\left(t,p\right)=x_0\left(t\right)+{{\displaystyle \sum }}_{k=1}^{+\infty }\hspace{0.17em}{x}_k\left(t\right)p^k$

(5)converges when $p=1$ . Therefore, the series solution of the form

$x\left(t\right)=x_0\left(t\right)+{{\displaystyle \sum }}_{n=1}^{+\infty }\hspace{0.17em}{x}_n\left(t\right),$

(6)which must satisfy the equation $N\left[x\left(t\right)\right]=0$ .

If we differentiate the zero-th order homotopy equation (3) by $m$ -times with respect to $p$ and dividing it by $m!$ finally if we set $p=0$ , we have the so-called $m^{\text{th}}$ order deformation equation as

$\mathcal{L}\left[x_m\left(t\right)-{\chi }_m{x}_{m-1}\left(t\right)\right]=R_{m-1}\left({\overrightarrow{x}}_{m-1}\left(t\right)\right)$

(7)where

$R_k\left(\overrightarrow{x}\left(t\right)\right)={\frac{1}{k!}\frac{{\partial }^k}{\partial p^k}\mathcal{N}\left[{{\displaystyle \sum }}_{m=1}^{+\infty }\hspace{0.17em}{x}_m\left(t\right)p^n\right]\hspace{0.17em}|}_{p=0}$

(8)and

${\chi }_k=\{\begin{array}{ll}0,& \hspace{0.17em}k\le 1;\hspace{0.17em}\\ 1,& \hspace{0.17em}k\hspace{0.17em}>1.\hspace{0.17em}\end{array}$

(9)Therefore, from (7) we have particular solution as

$x_m^{*}\left(t\right)={\chi }_m{x}_{m-1}\left(t\right)-{\mathcal{L}}^{-1}{R}_{m-1}\left({\overrightarrow{x}}_{m-1}\left(t\right)\right)$

(10)

Hence, the general solution becomes

$x_m\left(t\right)=x_m^{*}\left(t\right)+{\displaystyle \sum }_{k=1}^n\hspace{0.17em}{c}_k{f}_k\left(t\right)$

(11)where, $n$ is the order of given equation (1), $c_k$ are constant coefficients, and $f_k\left(t\right)$ are non-zero solution of $\mathcal{L}\left[x\left(t\right)\right]=0$ . As noted in form of the auxiliary linear operator equation (4), the linear operator $\mathcal{L}\left[x\left(t\right)\right]$ contains the auxiliary roots ${\lambda }_i$ of $\mathcal{L}\left[x\left(t\right)\right]=0$ , $i=1,2,..,n$ . Now, if all are not simultaneously equal to zero and known. Then, the initial approximation as well as finial series solution contains those unknowns ${\lambda }_i$ . To compute these unknown roots ${\lambda }_i$ , $i=1,2,..,n$ we minimize the square residual error $△\left({\lambda }_i\right)$ define as

$△\left({\lambda }_i\right)=\underset{\mathrm{\Omega }}{{\displaystyle \int }}\mathcal{\hspace{0.17em}}\mathcal{N}{\left[\hspace{0.17em}{{\displaystyle \sum }}_{k=0}^m\hspace{0.17em}{x}_k\left(t\right)\hspace{0.17em}\right]}^{2}dx$

(12)where $x_k\left(t\right)$ is the kth order homotopy approximation.

In order to reduce the computational time, we utilize the discretised form of the square residual which is based on newton quotes quadrature rule (simpson’s $\frac{1}{3}$ composite rule) as

$△\left({\lambda }_{\mathrm{i}}\right)\approx \frac{h}{3}{{\displaystyle \sum }}_{j=1}^{\frac{k}{2}}\hspace{0.17em}\left\{\mathcal{N}{\left[{{\displaystyle \sum }}_{i=0}^m\hspace{0.17em}{x}_i\left(t_{2j-2}\right)\right]}^2+4\mathcal{N}{\left[{{\displaystyle \sum }}_{i=0}^m\hspace{0.17em}{x}_i\left(t_{2j-1}\right)\right]}^2+\mathcal{N}{\left[{{\displaystyle \sum }}_{i=0}^m\hspace{0.17em}{x}_i\left(t_{2j}\right)\right]}^2\right\}$

(13)Where $f_k\left(\eta \right)$ is the homotopy approximation of the kth order, $h=\frac{b-a}{k}$ , $k=$ even number. Therefore, minimizeing the square residual $△\left({\lambda }_{\mathrm{i}}\right)$ , defined in equation (13), we compute unknowns ${\lambda }_{\mathrm{i}}$ . Thereby we ensure the convergence of our mth order series solution $x_m\left(t\right)$ (equation (11)). Readers are referred to our previous study ³³ for more detailed discussion on the necessary and sufficient condition for the convergence of the computed series solution.

Analytical and numerical illustration

Duffing oscillator

$\{\begin{array}{l}{\mathrm{x}}^{″}+\mathrm{x}+{\mathrm{x}}^3=0\\ \hspace{0.17em}\mathrm{x}\left(0\right)=1,{\mathrm{x}}^{\prime }\left(0\right)=0\end{array}$

(14)

As equation (14) is of 2^nd order differential equation, so from equation (4) we have the auxiliary linear operator of the form as

$\mathcal{L}\left[\mathrm{x}\left(\mathrm{t}\right)\right]=x^{″}-\left({\lambda }_1+{\lambda }_2\right)x^{\prime }+{\lambda }_1{\lambda }_{2}x$

(15)

For optimal homotopy perturbation based analytical solution at least one root ${\lambda }_1$ or ${\lambda }_2$ is to be remain unknown, which will be obtained optimally using residual minimization. Let ${\lambda }_1\ne 0$ be unknown and ${\lambda }_2=0\hspace{0.17em}.$ There are other numerious possibilities for the selection of auxiliary linear operator. One may keep both ${\lambda }_1$ and ${\lambda }_2$ as unknown. In that case computational time will be more.^31,32 Therefore, equation (15) becomes

$\mathcal{L}\left[\mathrm{x}\left(\mathrm{t}\right)\right]=x^{″}-{\lambda }_1{x}^{\prime }$

(16)Now, for the periodic solution we consider the initial approximation as $x_0\left(t\right)=A\hspace{0.17em}\mathit{cos}\left(\omega t\right)$ Using the mth order deformation equation defined in equation (7) with

$R_{m-1}\left({\overrightarrow{x}}_{m-1}\left(t\right)\right)={\mathrm{x}}_{\mathrm{m}-1}^{″}+x_{m-1}+{\displaystyle \sum }_n^{m-1}{x}_{m-n-1}{\displaystyle \sum }_{k=0}^n{x}_k\hspace{0.17em}{x}_{n-k}$

(17)and with the help of symbolic software we compute successive homotopy perturbation approximations. The first order approximation is

${\mathrm{x}}_1\left(\mathrm{t}\right)=\mathrm{A}\left(3{\mathrm{A}}^2\omega \hspace{0.17em}\cos \left[3\mathrm{t}\omega \right]\left({\mathrm{\omega }}^2+{\mathrm{\lambda }}_1^2\right)+{\mathrm{A}}^2\hspace{0.17em}\sin \left[3\mathrm{t}\omega \right]{\mathrm{\lambda }}_1\left({\mathrm{\omega }}^2+{\mathrm{\lambda }}_1^2\right)+3\mathrm{\omega }\left(4+3{\mathrm{A}}^2-4{\mathrm{\omega }}^2\right)\cos \left[\mathrm{t}\omega \right]\left(9{\mathrm{\omega }}^2+{\mathrm{\lambda }}_1^2\right)+3\left(4+3{\mathrm{A}}^2-4{\mathrm{\omega }}^2\right)\sin \left[\mathrm{t}\omega \right]{\mathrm{\lambda }}_1\left(9{\mathrm{\omega }}^2+{\mathrm{\lambda }}_1^2\right)-12{\mathrm{e}}^{\mathrm{t}{\mathrm{\lambda }}_1}\mathrm{\omega }\left({\mathrm{\omega }}^2\left(9+7{\mathrm{A}}^2-9{\mathrm{\omega }}^2\right)+\left(1+{\mathrm{A}}^2-{\mathrm{\omega }}^2\right){\mathrm{\lambda }}_1^2\right)\right)/12\mathrm{\omega }\left({\mathrm{\omega }}^2+{\mathrm{\lambda }}_1^2\right)\left(9{\mathrm{\omega }}^2+{\mathrm{\lambda }}_1^2\right)$

(18)Now, the condition for which it becomes periodic is that the coefficient of cos[ $\omega t]$ should be zero. Therefore,

$\{\hspace{0.17em}\begin{array}{l}\frac{\mathrm{A}\left(4+3{\mathrm{A}}^2-4{\mathrm{\omega }}^2\right)}{4\left({\mathrm{\omega }}^2+{\mathrm{\lambda }}_1^2\right)}=0\\ \frac{A\left(4+3A^2-4{\mathrm{\omega }}^2\right){\lambda }_1}{4\mathrm{\omega }\left({\mathrm{\omega }}^2+{\lambda }_1^2\right)}=0\end{array}$

(19)

Solving equation (19) we have amplitude $A=1$ and frequency $\mathrm{\omega }=\frac{1}{2}\sqrt{4+3{\mathrm{A}}^2}$ . Substituiting the value of amplitude A we have $\mathrm{\omega }=1.322875$ .

Now substituting the value of amplitude A and frequency $\mathrm{\omega }$ in the equation (20) we have the first order OM-HPM series solution as

$\mathrm{x}\left(\mathrm{t}\right)={\mathrm{x}}_0\left(\mathrm{t}\right)+x_1\left(t\right)$

(20)

$\mathrm{x}\left(\mathrm{t}\right)=\frac{3\sqrt{7}\left(-e^{t{\lambda }_1}+63\hspace{0.17em}\cos \left[\frac{\sqrt{7}t}{2}\right]+\mathrm{Cos}\left[\frac{3\sqrt{7}t}{2}\right]\right)+2\hspace{0.17em}\sin \left[\frac{3\sqrt{7}t}{2}\right]{\lambda }_1+12\sqrt{7}\cos \left[\frac{\sqrt{7}t}{2}\right]{\lambda }_1^2}{3\sqrt{7}\left(63+4{\lambda }_1^2\right)}$

(21)Therefore, using the minimization technique on $△\left({\lambda }_{\mathrm{i}}\right)$ in (13), we have the minimum square residual as $3.79662\times {10}^{-4}$ at ${\lambda }_1=-0.3078049$ . Hence the final solution becomes

$\mathrm{x}\left(\mathrm{t}\right)=-0.015778e^{-0.3078049t}+0.9999\hspace{0.17em}\cos \left[\frac{\sqrt{7}t}{2}\right]+0.0157781\hspace{0.17em}\cos \left[\frac{3\sqrt{7}t}{2}\right]-\hspace{0.17em}0.0012237\hspace{0.17em}\sin \left[\frac{3\sqrt{7}t}{2}\right]$

(22)

As the first order solution give us a sufficiently enough small residual error, the higher order may be avoided here. It is important to emphasized that all the analyses given in six have been considered as dimensionless. To confirm the solution accuracy we present Fig-1 where we compare the computed OM-HPM series based analytical solution with those of numerical solution obtained by 4^th order Rungee-Kutta method.

From Figure 1 it is clear that first order OM-HPM ³⁶ analytical series solution is almost matching with the numerical solution obtained by RK4 method. Therefore, for this problem first order OM-HPM approximation is sufficient. Then we continue upto 5^th order OM-HPM approximations. Table 1 represnts the different order of approximations, optimal value of ${\lambda }_1$ and the minimum square residual. From Table 1 one can easily see that the error diminishes with the order of approximations. Figure 2 is presented to show the phase portrait for velocity $x^{\prime }\left(t\right)$ versus the solution x(t) at different order of OM-HPM approximations. This problem was solved by He ³⁴ in 1998 to certify his proposed HPM and He obtained the relative percentage error upto 5.4%. It may be noted that there is no guarantee of convergence of the solution by HPM .^{3,8,13–20,27,29,32} Then Tao et al. ⁹ applied VIM with matrix Lagrange’s multiplier to solve this duffing oscillator and obtained a series based semi-analytical solution but there was no improvement on the solution. While by our recently proposed OM-HPM ³⁶ one can get accurate and convergent solution, and this efficiency of OM-HPM and inefficiency of HPM/VIM lies on the choice of optimal linear operator. We found the maximum relative percentage error as 1.303%. We also presented the square residual in Table 1 to clarify the high efficiency of our advanced technique OM-HPM ³⁶ as order increases.

Figure 1.

Comparison of the computed 1^st order OM-HPM ³² analytical solution with the numerical solution by RK4. (a). Time series solution of x(t) versus t. (b). Time series of velocity $x^{\prime }\left(t\right)$ versus t.

Table 1.

Optimal table for duffing oscillator (14).

Order m	Optimal ${\lambda }_1$	Min error
1	$-0.3078049$	$3.79662\times {10}^{-4}$
3	$-0.56026707$	$8.54681\times {10}^{-6}$
5	$-1.478902$	$1.23794\times {10}^{-7}$

Figure 2.

Phase portrait for velocity $x^{\prime }\left(t\right)$ versus the solution x(t) for the equation (14)) at different order of approximations.

From Figure 2 it is clear that both the lower (m = 1) and higher (m = 5) order OM-HPM approximation are close to each other and providing limit cycle, that is, stable periodic motion over time, forming a closed trajectory in the phase space.

Van der Pol oscillator

$\{\begin{array}{l}{x}^{″}\left(t\right)+x\left(t\right)-x^{\prime }\left(t\right)\left\{1-b{x}^2\left(t\right)\right\}=0\\ \mathrm{x}\left(0\right)=1,{\mathrm{x}}^{\prime }\left(0\right)=0\end{array}$

(23)

Several researchers attempted to solve this van der pol oscillator analytically using the HPM,^13,34,35 HAM,^40,41 and VIM.^11,35 Ozis et al. ³⁴ and others ^13,35 solved using the classical HPM and obtained series solution with relative percentage error upto 37.5%. Here, we are interested to apply our recently proposed advanced technique OM-HPM. ³⁴ Similarly, using the linear operator (17), the initial approximation $x_0\left(t\right)=A\hspace{0.17em}\mathit{cos}\left(\omega t\right),$ and the mth order deformation equation as

$R_{m-1}\left({\overrightarrow{x}}_{m-1}\left(t\right)\right)={\mathrm{x}}_{\mathrm{m}-1}^{″}+x_{m-1}-x_{m-1}^{\prime }+b{{\displaystyle \sum }}_n^{m-1}{x}_{m-n-1}^{\prime }{{\displaystyle \sum }}_{k=0}^n{x}_k\hspace{0.17em}{x}_{n-k}$

(24)

we compute successive homotopy perturbation approximations. The first order approximation is of the form

${\mathrm{x}}_1\left(\mathrm{t}\right)=\frac{1}{12\omega {\lambda }_1\left({\omega }^2+{\lambda }_1^2\right)\left(9{\omega }^2+{\lambda }_1^2\right)}A(36\left(-3+A^{2}b\right)\left(-1+e^{t{\lambda }_1}\right){\omega }^5-3{\omega }^3\left(36e^{t{\lambda }_1}-36e^{t{\lambda }_1}{\omega }^2+36\left(-1+{\omega }^2\right)\cos \left[t\omega \right]+9(-4+A^{2}b\right)\omega \hspace{0.17em}\sin \left[t\omega \right]+A^{2}b\omega \hspace{0.17em}\sin \left[3t\omega \right]){\lambda }_1-{\omega }^2(-27(-4+A^{2}b)\omega \hspace{0.17em}\cos \left[t\omega \right]-A^{2}b\omega \hspace{0.17em}\cos \left[3t\omega \right]+4(-(30-3e^{t{\lambda }_1}+A^{2}b\left(-10+3e^{t{\lambda }_1}\right))\omega +27\left(-1+{\omega }^2\right)\sin [t\omega ])){\lambda }_1^2-3\omega \left(4e^{t{\lambda }_1}-4e^{t{\lambda }_1}{\omega }^2+4\left(-1+{\omega }^2\right)\cos \left[t\omega \right]+(-4+A^{2}b\right)\omega \hspace{0.17em}\sin \left[t\omega \right]+A^{2}b\omega \hspace{0.17em}\sin [3t\omega ]){\lambda }_1^3+\left(3\left(-4+A^{2}b\right)\omega \hspace{0.17em}\cos \left[t\omega \right]+A^{2}b\omega \hspace{0.17em}\cos \left[3t\omega \right]-4((-3+A^{2}b\right)\omega +3(-1+{\omega }^2)\sin [t\omega ])){\lambda }_1^4$

(25)

Therefore, for the periodic solution we have

$\{\begin{array}{l}\frac{A\left(4-4{\omega }^2+\left(-4+A^{2}b\right){\lambda }_1\right)}{4\left({\omega }^2+{\lambda }_1^2\right)}=0\\ -\frac{A\left(\left(-4+A^{2}b\right){\omega }^2+4\left(-1+{\omega }^2\right){\lambda }_1\right)}{4\omega \left({\omega }^2+{\lambda }_1^2\right)}=0\end{array}$

(26)Solving we get amplitude $A=\frac{2}{\sqrt{b}}$ and frequency $\omega =1$ . For b = 1 the limiting amplitude is A = 2 and frequency $\omega$ as 1. Substituiting the value of computed amplitude A and frequency $\omega ,$ we have the first order OM-HPM series solution as

$\mathrm{x}\left(\mathrm{t}\right)=\frac{2\left(9\left(-1+e^{t{\lambda }_1}\right)-3\left(-9\hspace{0.17em}\cos \left[t\right]+\sin \left[3t\right]\right){\lambda }_1+\left(-1+\cos \left[3t\right]\right){\lambda }_1^2+3\hspace{0.17em}\cos \left[t\right]{\lambda }_1^3\right)}{3{\lambda }_1\left(9+{\lambda }_1^2\right)}$

(27)Therefore, we have the minimum square residual for b = 1 as $4.066553\times {10}^{-2}$ at the ${\lambda }_1=-2.4889329.$ Again substituting the optimal value of the auxiliary root ${\lambda }_1$ into the equation (27) we have the final solution as

$\mathrm{x}\left(\mathrm{t}\right)=-0.009989\left(9\left(-1+e^{-3.3262438t}\right)-110.4036\mathrm{Cos}\left[t\right]+11.063897(-1+\cos \left[3t\right])+9.97873(-9\hspace{0.17em}\cos \left[t\right]+\sin \left[3t\right]\right))$

(28)

However, here the residual error at 1^st order is not satisfactory. Here, we continue upto 5^th order approximations and prepare the optimal table as

Therefore, from Table 2 we have the best optimal auxiliary linear operator as.

$\mathcal{L}\left[\mathrm{x}\left(\mathrm{t}\right)\right]=x^{″}+1.2530439x^{\prime }$

(29)

Table 2.

Optimal table for Van der Pol oscillator (23).

Order m	Optimal ${\lambda }_1$	Min error
1	$-2.4889329$	$4.066553\times {10}^{-2}$
3	$-1.8649311$	$7.37648\times {10}^{-4}$
5	$-1.2530439$	$2.06427\times {10}^{-5}$

Now, substituting the optimal value of the auxiliary roots ${\lambda }_1$ we have the 5^th order approximation solution of the form

$\mathrm{x}\left(\mathrm{t}\right)={\mathrm{x}}_0\left(\mathrm{t}\right)+x_1\left(t\right)+x_2\left(t\right)+{\mathrm{x}}_3\left(\mathrm{t}\right)+x_4\left(t\right)+x_5\left(t\right)$

(30)

The solutions of the equation (23) are shown in the Figure 3. In the Figure 3 we compare our computed OM-HPM solution with the numerical solution by RK4. Figure 3(a) presents the solution x(t) while Figure 3(b) represents the velocity $x^{\prime }\left(t\right)$ . From Figure 3(a) to (b) it is seen that our 1^st order approximation solution is not much well but our 3^rd and 5^th order approximations are very close to the numerical solution and we have the relative error as 3.04%, 1.93% and 0.84% for 1^st, 3^rd and 5^th approximations, respectively. Eigoli et al ⁴¹ solved this oscillator using HAM and found that to achieve good accuracy one need 25^th order approximations. Moreover, when Huang ³⁵ applied VIM for this oscillator, he found that VIM gives divergent solution. Therefore, Huang ³⁵ improved the scheme for the VIM by adding HPM^3,13,14 and Laplace transformation and shown that 6^th order solution by improved VIM is agree with numerical solution. On the other hand, from Table 2 and Figure 3(a) to (b) it is noteworthy that our OM-HPM ³⁴ gives better solution at only 1^st order approximation. Therefore, in comparison to the classical HPM, HAM, ⁴¹ VIM, ³⁵ and improved VIM, ³⁵ our OM-HPM ³⁶ is much more efficient to handle such strongly nonlinear oscillators.

Figure 3.

Comparison of the computed analytical solution with the numerical solution by RK4 at different order of approximations (m = 1,3 and 5) for the limiting case b = 1. (a). Time series solution of x(t) (b). Time series of velocity $x^{\prime }\left(t\right)$ .

The phase portrait at the different order (m = 1 and 5) of OM-HPM approximations series solutions for the oscillator (23) is shown in the Figure 4. We have seen that there is a closed trajectory surrounding the origin. Hence, the system exhibits limit cycle behavior. The system does not settle down to equilibrium but rather to a periodic solution.

Figure 4.

Phase portrait for velocity $x^{\prime }\left(t\right)$ versus the solution x(t) for the equation (23) at different order of approximations.

Rayleigh oscillator

$\{\begin{array}{l}{x}^{″}+x-x^{\prime }\hspace{0.17em}\left(1-{x^{\prime }}^2\right)=0\\ x\left(0\right)=1,x^{\prime }\hspace{0.17em}\left(0\right)=0\end{array}$

(31)

Following the similar process with the linear operator (17), the initial approximation $x_0\left(t\right)=A\hspace{0.17em}\mathit{cos}\left(\omega t\right)$ and the mth order deformation equation

$R_{m-1}\left({\overrightarrow{x}}_{m-1}\left(t\right)\right)={\mathrm{x}}_{\mathrm{m}-1}^{″}+x_{m-1}-x_{m-1}^{\prime }+{\displaystyle \sum }_n^{m-1}{x}_{m-n-1}^{\prime }{\displaystyle \sum }_{k=0}^n{x}_k^{\prime }\hspace{0.17em}{x}_{n-k}^{\prime }$

(32)we compute successive homotopy perturbation approximations. The first order approximation is

${\mathrm{x}}_1\left(\mathrm{t}\right)=\frac{1}{12\omega {\lambda }_1\left({\omega }^2+{\lambda }_1^2\right)\left(9{\omega }^2+{\lambda }_1^2\right)}A(36\left(-1+e^{t{\lambda }_1}\right){\omega }^5\left(-3+2A^2{\omega }^2\right)+3{\omega }^3\left(-36e^{t{\lambda }_1}+36e^{t{\lambda }_1}{\omega }^2-36\left(-1+{\omega }^2\right)\cos \left[t\omega \right]-9\omega (-4+3A^2{\omega }^2\right)\sin \left[t\omega \right]+A^2{\omega }^3\hspace{0.17em}\sin [3t\omega ]){\lambda }_1-{\omega }^2(-27\omega \left(-4+3A^2{\omega }^2\right)\cos \left[t\omega \right]+A^2{\omega }^3\hspace{0.17em}\cos \left[3t\omega \right]+4(\omega (-30+3e^{t{\lambda }_1}+20A^2{\omega }^2)+27\left(-1+{\omega }^2\right)\sin [t\omega ])){\lambda }_1^2+3\omega \left(-4e^{t{\lambda }_1}+4e^{t{\lambda }_1}{\omega }^2-4\left(-1+{\omega }^2\right)\cos \left[t\omega \right]+(4\omega -3A^2{\omega }^3\right)\sin \left[t\omega \right]+A^2{\omega }^3\hspace{0.17em}\sin [3t\omega ]){\lambda }_1^3+\left(3\omega \left(-4+3A^2{\omega }^2\right)\cos \left[t\omega \right]-A^2{\omega }^3\hspace{0.17em}\cos \left[3t\omega \right]+4(3\omega -2A^2{\omega }^3-3\left(-1+{\omega }^2\right)\sin [t\omega ]\right)){\lambda }_1^4)$

(33)

Similarly for the periodic solution we have

$\{\begin{array}{l}\frac{A\left(4-4{\omega }^2+\left(-4+3A^2{\omega }^2\right){\lambda }_1\right)}{4\left({\omega }^2+{\lambda }_1^2\right)}=0\\ \frac{A\left(4{\omega }^2-3A^2{\omega }^4-4\left(-1+{\omega }^2\right){\lambda }_1\right)}{4\omega \left({\omega }^2+{\lambda }_1^2\right)}=0\end{array}$

(34)

Solving we have

$A=\frac{2}{\sqrt{3}}\hspace{0.17em}\text{and}\hspace{0.17em}\omega =1$

Substituting the value of computed amplitude A and frequency $\omega$ we have the first order OM-HPM series solution as

$\mathrm{x}\left(\mathrm{t}\right)={\mathrm{x}}_0\left(\mathrm{t}\right)+x_1\left(t\right)$

(35)

$\mathrm{x}\left(\mathrm{t}\right)=\frac{2\hspace{0.17em}\cos \left[t\right]}{\sqrt{3}}+\frac{-12\left(-1+e^{t{\lambda }_1}\right)+4\hspace{0.17em}\sin \left[3t\right]{\lambda }_1-\left(4\left(-\frac{10}{3}+3e^{t{\lambda }_1}\right)+\frac{4}{3}\cos \left[3t\right]\right){\lambda }_1^2+4\hspace{0.17em}\sin \left[3t\right]{\lambda }_1^3+\left(\frac{4}{3}-\frac{4}{3}\cos \left[3t\right]\right){\lambda }_1^4}{6\sqrt{3}{\lambda }_1\left(1+{\lambda }_1^2\right)\left(9+{\lambda }_1^2\right)}$

(36)

The minimum square residual is computed as $1.1618723\times {10}^{-3}$ at the ${\lambda }_1=-1.77395597.$ Again substituiting the optimal value of the auxiliary root ${\lambda }_1$ into the equation (36) we have final solution as

$\mathrm{x}\left(\mathrm{t}\right)-0.07232426+0.053587e^{-1.77395597t}+1.1547005\hspace{0.17em}\cos \left[t\right]+0.0187371499\hspace{0.17em}\cos \left[3t\right]+0.03168706\hspace{0.17em}\sin \left[3t\right]$

(37)In similar process we successively compute up to 5^th order OM-HPM series solutions and the solutions are shown in Figures 5 and 6. In Figure 5(a) we compare our computed OM-HPM position solution with numerical RK4 solution at different order of approximations while in Figure 5(b) we compare velocity. The phase portrait of the oscillator (31) is shown in Figure 6 at the different order of approximation solutions.

Figure 5.

Comparison of the computed analytical solution of equation (33) with the numerical solution by RK4 at different order of approximations (a) time series solution of x(t) (b) time series of velocity $x^{\prime }\left(t\right)$ .

Figure 6.

Phase portrait for velocity $x^{\prime }\left(t\right)$ versus the solution x(t) for the equation (24) at different order of approximations.

Similar to the Van der Pol oscillator, here also arise a closed trajectory surrounding the origin. Hence, the system exhibits limit cycle, that is, stable periodic motion over time.

Conclusions

In this paper, we applied the optimal and modified homotopy perturbation method (OM-HPM) technique to provide convergent semi-analytical solutions for three different models: duffing oscillator (DO), Rayleigh oscillator (RO) and van der Pol oscillator (vdPO). The high accuracy of this technique is attributed to the selection of optimal auxiliary linear operator, which makes it a practical tool for engineers and researchers. It has been shown that out of the three models discussed, only the DO is conservative in nature and other two are non-conservative in nature. Surprisingly all the three oscillators are stable oscillators. Our method is suitable as long as the phases remains closed in nature, thereby justifying the stability. Present method will not give encouraging results on unstable cases reflecting open nature of phase portrait. It may be recommended that this advanced approach OM-HPM can be applied to other strongly nonlinear oscillators to explore the complex dynamics. However, OM-HPM need more CPU time than traditional HPM.