A simple modified harmonic balance method for strongly nonlinear oscillator with cubic non-linearity and harmonic restoring force

Abstract

In this article, a very simple modified form of the harmonic balance method is used to solve a strongly nonlinear oscillator with cubic nonlinearity and harmonic restoring force. Taylor series expansion up to third term is considered for the harmonic restoring force. The first approximate solutions of the present method pleasantly agree with the numerical solution obtained by Runge–Kutta fourth order method. Accuracy and simplicity of the present method solution is established when compared with the other method solutions. The present method can be utilized to other nonlinear oscillators.

Introduction

Nonlinear oscillators are formed from different aspects of engineering problems. Nonlinear oscillators can be modeled by ordinary and partial differential equations. Solution of the nonlinear ordinary and partial differential equations are not easy like linear differential equations. Moreover, exact solutions of the nonlinear oscillators are extremely complex and most of the times unreachable. Therefore, many researchers have contributed to obtain the approximate solution of various nonlinear oscillators. Initially the perturbation methods^1–3 were used by expanding a small parameter to solve nonlinear oscillatory systems. But such expansion was not easy and suitable. Latter homotopy perturbation method ⁴ was first proposed by Ji-Huan He. Many approximate methods such as the harmonic balance method, ⁵ the energy balance method, ⁶ He’s frequency formulation method, ⁷ the variational iteration method, ⁸ the iteration perturbation method, ⁹ the homotopy analysis method, ¹⁰ max-min approach, ¹¹ VIM-Pade technique ¹² etc. were proposed by many researchers to solve strongly nonlinear oscillators.

Homotopy perturbation method was applied by others^13,14 appropriately. Also, many researchers were applied the harmonic balance method, ¹⁵ the energy balance method, ¹⁶ He’s frequency formulation method,^17,18 the variational iteration method,^19–22 the homotopy analysis method,^23,24 VIM-Pade technique^25,26 etc. successfully to solve strongly nonlinear problems.

Many types of modification like the global residue harmonic balance method,^27–31 the modified energy balance method,^32,33 the modified harmonic balance method,^34–36 the modified homotopy perturbation method,^37–40 the modified Lindested–Poincare method,^41,42 the modified Laplace based variational iteration method ⁴³ etc. were also used by many researchers to solve strong nonlinear problems. Coupling techniques between the homotopy perturbation method and the variational method ⁴⁴ and between the Lindested–Poincare method and the homotopy perturbation method ⁴⁵ were also used for solving the nonlinear oscillators. Harmonic balance method and different modified forms of the harmonic balance method have been widely used for solving nonlinear oscillators from last few decades. In harmonic balance method the response of the system is assumed in terms of a Fourier series and using this expression in the governing nonlinear differential equation and separating the coefficients of the harmonic terms one can obtain the unknown coefficients and frequency amplitude relation of the nonlinear system. When the HBM is applied to the nonlinear equations for higher-order approximation, then a set of difficult nonlinear complex equations appear and it is very difficult to analytically solve these complex equations. In the present modified form of harmonic balance method, we applied a simple factor on the strongly nonlinear oscillator. To find the higher-order approximate solution and the corresponding frequencies, Fourier series expansion is used and the set of algebraic equations arises is solved easily.

In the present article, we have used the simple modified form of the harmonic balance method ³⁶ to solve a strongly nonlinear oscillator with cubic nonlinearity and harmonic restoring force modeled by ³⁰

$\ddot{x}+x+a{x}^3+b\sin x=0,$

(1)

where a and b are known constants and initial conditions are $x\left(0\right)=A,\dot{x}\left(0\right)=0$ .

The nonlinear equation represented by equation (1) has application in mechanical engineering (see ⁴⁶ for details) and transmission system consisting of porous media. This strongly nonlinear oscillator can be modeled for the vibration of the system consisting a mass resting on a spring with cubic nonlinearity and the driving force as illustrated by Figure 1. Here we consider the system response, mass and the linear spring stiffness coefficient by $x\left(t\right),M$ and $k$ , respectively. Also, the term $\sin x$ resulting for the driving force.

Figure 1.

Graphical configuration of the problem.

Marinca and Herisanu ⁴⁷ applied an iteration method to solve the oscillator represented by equation (1). Also, many analytical techniques have been used to solve the strongly nonlinear oscillator with cubic nonlinearity and harmonic restoring force like Lu and Ma ²⁶ used the VIM-pade technique, Alal et al. ³² proposed the modified energy balance method, Ismail and Alal ²⁹ applied the global residue harmonic balance method, Alal and Chowdhury ¹⁵ used the harmonic balance method and Junfeng Lu ³⁰ used the global residue harmonic balance method since the last decade.

The solution process for the equation (1) using the harmonic balance method and the energy balance method is laborious and complex. Also, the solution deteriorates for larger amplitudes and time. Solution obtained by the global residue harmonic balance method needed second or larger approximations to improve the result. Therefore, a simple and efficient harmonic balance method ³⁶ presented in this article to overcome the above-mentioned problems. Only the first approximate solution of the present method provides better results than the other existing methods.

Modified harmonic balance method

Let us consider a strongly nonlinear oscillator

$\ddot{x}+{\omega }_0^2\hspace{0.17em}x+f\left(x\right)=0.$

(2)

where $f\left(x\right)$ is a nonlinear function such that $f\left(-x\right)=-f\left(x\right)$ , ${\omega }_0\ge 0$ and over dots represents derivative with respect to $t$ . The initial conditions are given as $x\left(0\right)=A,\dot{x}\left(0\right)=0$ .

The n^th order approximate solution of equation (2) is considered as ³¹

$x\left(t\right)=A\left(\left(1-u-v-\dots \right)\cos \phi +u\hspace{0.17em}\cos 3\hspace{0.17em}\phi +v\hspace{0.17em}\cos 5\hspace{0.17em}\phi +\dots \right),$

(3)

where $u,v,\cdots$ are unknown parameters to be determined later and $\phi =\omega \hspace{0.17em}t$ , $\omega$ is an unknown angular frequency.

The zeroth order approximate solution of the equation (2) is considered choosing $u=v=\cdots =0$

$x_0\left(t\right)=A\hspace{0.17em}\cos \phi .$

(4)

Using equation (4) in equation (2) and using the condition of no secular term, $\omega$ can be found.

According to the present modified form of harmonic balance method, we rewrite the equation (2) to obtain the higher order solutions as

$\frac{{\omega }^2{x}^{″}+{\omega }_0^2\hspace{0.17em}x+f\left(x\right)}{1+{x_0}^2}=0.$

(5)

where prime denotes differentiation with respect to $\phi$ . It is worth mentioning that dividing the equation (2) by the factor $\left(1+{x_0}^2\right)$ and using the present modified form of the harmonic balance method, a set of simple algebraic equations appear which contains lower order terms and these lower order terms form the solution promptly convergent. On the other hand, dividing by the factor $\left(1+{x_0}^2\right)$ , the equation (2) does not lose its originality.

The initial conditions are replace by $x\left(0\right)=A,x^{\prime }\left(0\right)=0$ .

Then equation (5) is expanded in Fourier series

$\frac{{\omega }^2{x}^{″}+{\omega }_0^2\hspace{0.17em}x+f\left(x\right)}{1+{x_0}^2}=m_1\left(\omega ,u,v,\cdots \right)\cos \phi +m_3\left(\omega ,u,v,\cdots \right)\cos 3\phi +m_5\left(\omega ,u,v,\cdots \right)\cos 5\phi +\cdots \cdots \cdots ,$

(6)

where

$m_{2n-1}\left(\omega ,u,v,\cdots \right)=\frac{4}{\pi }{\int }_0^{\frac{\pi }{2}}\left(\frac{{\omega }^2{x}^{″}+{\omega }_0^2\hspace{0.17em}x+f\left(x\right)}{1+{x_0}^2}\right)\cos \left(2n-1\right)\phi \hspace{0.17em}d\phi ,n=0,1,2,\cdots ,$

(7)

after putting the values of $x$ and $x_0$ given by Equations (3) and (4), respectively.

From Equations (5) and (6), we get

$m_1\left(\omega ,u,v,\cdots \right)\cos \phi +m_3\left(\omega ,u,v,\cdots \right)\cos 3\phi +m_5\left(\omega ,u,v,\cdots \right)\cos 5\phi +\dots =0.$

(8)

Now equating the coefficients of $\cos \phi$ , $\cos 3\phi ,\cos 5\phi ,\dots$ from both sides of equation (8) we get a set of nonlinear algebraic equations involving the unknown parameters $u,v,\cdots$ along with unknown frequency $\omega$ . Neglecting the higher order terms of $u,v,\cdots$ and solving these algebraic equations, we can find out the values of $u,v,\cdots$ along with the frequency $\omega$ .

For the first approximate solution, we consider

$x\left(t\right)=A\left(\left(1-u\right)\cos \phi +u\hspace{0.17em}\cos 3\phi \right).$

(9)

Now substituting the values of $x$ and $x_0$ given by Equations (9) and (4), respectively, in the equation (5) and then expanding in a Fourier series and neglecting the higher harmonic terms more than third harmonic, we get

$\frac{{\omega }^2{x}^{″}+{\omega }_0^2\hspace{0.17em}x+f\left(x\right)}{1+{x_0}^2}=m_1\left(\omega ,u\right)\cos \phi +m_3\left(\omega ,u\right)\cos 3\phi .$

(10)

From Equations (5) and (10) we have

$m_1\left(\omega ,u\right)\cos \phi +m_3\left(\omega ,u\right)\cos 3\phi =0.$

(11)

Now equating the coefficient of $\cos \phi$ and $\cos 3\phi$ from both the sides of equation (11), we get

$m_1\left(\omega ,u\right)=0,$

(12)

$m_3\left(\omega ,u\right)=0.$

(13)

Neglecting higher order terms of $u$ and solving above two equations, the values of $u$ and $\omega$ can be found. Hence the first approximation solution is

$x\left(t\right)=A\left(\left(1-u\right)\cos \phi +u\hspace{0.17em}\cos 3\phi \right),$

(14)

where $u$ and $\omega$ can be found by solving Equations (12) and (13).

The solution

We consider the nonlinear oscillator

$\ddot{x}+x+a\hspace{0.17em}{x}^3+b\hspace{0.17em}\sin x=0.$

(15)

where ${\omega }_0^2=1$ and the initial conditions are given as $x\left(0\right)=A,\dot{x}\left(0\right)=0$ . From above equation, it is obvious that $f\left(-x\right)=-f\left(x\right)$ . Now one can expand $\sin x$ in a Taylor’s series expansion, $\sin x=x-\frac{x^3}{3!}+\frac{x^5}{5!}-\cdots$ . Neglecting the terms more than $x^5$ of $\sin x$ series equation (15) then reduces to

$\ddot{x}+x+a\hspace{0.17em}{x}^3+b\hspace{0.17em}\left(x-\frac{x^3}{3!}+\frac{x^5}{5!}\right)=0.$

(16)

Equation (16) can be written as

$\frac{{\omega }^2{x}^{″}+x+a\hspace{0.17em}{x}^3+b\hspace{0.17em}\left(x-\frac{x^3}{3!}+\frac{x^5}{5!}\right)}{1+{x_0}^2}=0.$

(17)

For the first approximate solution putting the values of $x_0\left(t\right)=A\hspace{0.17em}\cos \phi$ and $x\left(t\right)$ from equation (9) into equation (17) and then expanding in Fourier series and neglecting the higher harmonic terms, we get

$\frac{{\omega }^2{x}^{″}+x+a\hspace{0.17em}{x}^3+b\hspace{0.17em}\left(x-\frac{x^3}{3!}+\frac{x^5}{5!}\right)}{1+{x_0}^2}=m_1\left(\omega ,u\right)\cos \phi +m_3\left(\omega ,u\right)\cos 3\phi .$

(18)

Hence from Equations (17) and (18), we get

$m_1\left(\omega ,u\right)\cos \phi +m_3\left(\omega ,u\right)\cos 3\phi =0,$

(19)

where

$m_1\left(\omega ,u\right)={\mathrm{b}}_1+{\mathrm{b}}_2{\omega }^2+\left({\mathrm{b}}_3+{\mathrm{b}}_4{\omega }^2\right)u+O\left(u^2\right),$

$m_3\left(\omega ,u\right)={\mathrm{d}}_1+{\mathrm{d}}_2{\omega }^2+\left({\mathrm{d}}_3+{\mathrm{d}}_4{\omega }^2\right)u+d_5\hspace{0.17em}{u}^2+O\left(u^3\right),$

(20)

and

${\mathrm{b}}_1=\frac{2}{A}-\frac{2a}{A}+aA-\frac{2}{A\sqrt{1+A^2}}+\frac{2a}{A\sqrt{1+A^2}}+\frac{47b}{20A}-\frac{7Ab}{40}+\frac{A^{3}b}{160}-\frac{47b}{20A\sqrt{1+A^2}},$

${\mathrm{b}}_2=-\frac{2}{A}+\frac{2}{A\sqrt{1+A^2}},$

${\mathrm{b}}_3=-\frac{8}{A^3}+\frac{24a}{A^3}-\frac{4}{A}+\frac{12a}{A}-3aA+\frac{8}{A^3\sqrt{1+A^2}}-\frac{24a}{A^3\sqrt{1+A^2}}+\frac{8}{A\sqrt{1+A^2}}-\frac{24a}{A\sqrt{1+A^2}}-\frac{37b}{3A^3}-\frac{37b}{6A}+\frac{13Ab}{24}-\frac{A^{3}b}{48}+\frac{37b}{3A^3\sqrt{1+A^2}}+\frac{37b}{3A\sqrt{1+A^2}},$

${\mathrm{b}}_4=\frac{72}{A^3}+\frac{20}{A}-\frac{72}{A^3\sqrt{1+A^2}}-\frac{56}{A\sqrt{1+A^2}},$

${\mathrm{d}}_1=-\frac{8}{A^3}+\frac{8a}{A^3}-\frac{2}{A}+\frac{2a}{A}+\frac{8}{A^3\sqrt{1+A^2}}-\frac{8a}{A^3\sqrt{1+A^2}}+\frac{6}{A\sqrt{1+A^2}}-\frac{6a}{A\sqrt{1+A^2}}-\frac{47b}{5A^3}-\frac{47b}{20A}+\frac{A^{3}b}{480}+\frac{47b}{5A^3\sqrt{1+A^2}}+\frac{141b}{20A\sqrt{1+A^2}},$

${\mathrm{d}}_2=\frac{8}{A^3}+\frac{2}{A}-\frac{8}{A^3\sqrt{1+A^2}}-\frac{6}{A\sqrt{1+A^2}},$

${\mathrm{d}}_3=\frac{32}{A^5}-\frac{96a}{A^5}+\frac{40}{A^3}-\frac{120a}{A^3}+\frac{8}{A}-\frac{24a}{A}+3aA-\frac{32}{A^5\sqrt{1+A^2}}+\frac{96a}{A^5\sqrt{1+A^2}}-\frac{56}{A^3\sqrt{1+A^2}}+\frac{168a}{A^3\sqrt{1+A^2}}-\frac{24}{A\sqrt{1+A^2}}+\frac{72a}{A\sqrt{1+A^2}}+\frac{148b}{3A^5}+\frac{185b}{3A^3}+\frac{37b}{3A}-\frac{13Ab}{24}+\frac{A^{3}b}{96}-\frac{148b}{3A^5\sqrt{1+A^2}}-\frac{259b}{3A^3\sqrt{1+A^2}}-\frac{37b}{A\sqrt{1+A^2}},$

${\mathrm{d}}_4=-\frac{288}{A^5}-\frac{296}{A^3}-\frac{56}{A}+\frac{288}{A^5\sqrt{1+A^2}}+\frac{440}{A^3\sqrt{1+A^2}}+\frac{168}{A\sqrt{1+A^2}},$

${\mathrm{d}}_5=\frac{384a}{A^7}+\frac{864a}{A^5}+\frac{576a}{A^3}+\frac{84a}{A}-9aA-\frac{384a}{A^7\sqrt{1+A^2}}-\frac{1056a}{A^5\sqrt{1+A^2}}-\frac{960a}{A^3\sqrt{1+A^2}}-\frac{288a}{A\sqrt{1+A^2}}-\frac{224b}{3A^7}-\frac{168b}{A^5}-\frac{112b}{A^3}-\frac{49b}{3A}+\frac{7Ab}{4}-\frac{A^{3}b}{16}+\frac{224b}{3A^7\sqrt{1+A^2}}+\frac{616b}{3A^5\sqrt{1+A^2}}+\frac{560b}{3A^3\sqrt{1+A^2}}+\frac{56b}{A\sqrt{1+A^2}}.$

(21)

Now equating the coefficient of $\cos \phi$ and $\cos 3\phi$ from both sides of equation (19), and neglecting the higher order terms of $u$ , we get

${\mathrm{b}}_1+{\mathrm{b}}_2{\omega }^2+\left({\mathrm{b}}_3+{\mathrm{b}}_4{\omega }^2\right)u=0,$

(22)

${\mathrm{d}}_1+{\mathrm{d}}_2{\omega }^2+\left({\mathrm{d}}_3+{\mathrm{d}}_4{\omega }^2\right)u+{\mathrm{d}}_5\hspace{0.17em}{u}^2=0.$

(23)

Eliminating ${\omega }^2$ from Equations (22) and (23) and solving for $u$ , we get

$u=\frac{-{\mathrm{b}}_4{\mathrm{d}}_1+{\mathrm{b}}_3{\mathrm{d}}_2-{\mathrm{b}}_2{\mathrm{d}}_3+{\mathrm{b}}_1{\mathrm{d}}_4+\sqrt{{\left({\mathrm{b}}_4{\mathrm{d}}_1-{\mathrm{b}}_3{\mathrm{d}}_2+{\mathrm{b}}_2{\mathrm{d}}_3-{\mathrm{b}}_1{\mathrm{d}}_4\right)}^2-4\left({\mathrm{b}}_2{\mathrm{d}}_1-{\mathrm{b}}_1{\mathrm{d}}_2\right)\left({\mathrm{b}}_4{\mathrm{d}}_3-{\mathrm{b}}_3{\mathrm{d}}_4+{\mathrm{b}}_2{\mathrm{d}}_5\right)}}{2\left({\mathrm{b}}_4{\mathrm{d}}_3-{\mathrm{b}}_3{\mathrm{d}}_4+{\mathrm{b}}_2{\mathrm{d}}_5\right)}.$

(24)

Solving equation (22) for $\omega$ , we get

$\omega =\frac{\sqrt{-{\mathrm{b}}_1-{\mathrm{b}}_{3}u}}{\sqrt{{\mathrm{b}}_2+{\mathrm{b}}_{4}u}}.$

(25)

Hence the first approximation solution of equation (15) is

$x\left(t\right)=A\left(\left(1-u\right)\cos \phi +u\hspace{0.17em}\cos 3\phi \right),$

(26)

where $u$ and $\omega$ is given by equations (24) and (25), respectively.

He’s frequency formulation method solution

Considering two trial solutions ⁷ $x_1=A\cos {\omega }_{1}t$ and $x_2=A\cos {\omega }_{2}t$ for the equation (16) we have obtained the residuals as

$R_1\left(t\right)=\frac{1}{120}A{\mathit{cos}\omega }_{1}t\left(120\left(1+b-{{\mathrm{\omega }}_1}^2\right)+20A^2\left(6a-b\right){\cos }^2{\mathrm{\omega }}_{1}t+A^{4}b\hspace{0.17em}{\cos }^4{\mathrm{\omega }}_{1}t\right),$

(27)

$R_2=\frac{1}{120}A{\mathit{cos}\omega }_{2}t\left(120\left(1+b-{{\mathrm{\omega }}_2}^2\right)+20A^2\left(6a-b\right){\cos }^2{\mathrm{\omega }}_{2}t+A^{4}b\hspace{0.17em}{\cos }^4{\mathrm{\omega }}_{2}t\right).$

(28)

We have found the values of $\stackrel{\sim }{R_1}$ and $\stackrel{\sim }{R_2}$ as given below

$\stackrel{\sim }{R_1}=\frac{1}{T}{\int }_0^T{R}_1\left(t\right){\mathit{cos}\omega }_{1}tdt=\frac{1}{384}A\left(144a{A}^2+\left(192-24A^2+A^4\right)b-192\left(-1+{{\mathrm{\omega }}_1}^2\right)\right),$

(29)

$\stackrel{\sim }{R_2}=\frac{1}{T}{\int }_0^T{R}_1\left(t\right){\mathit{cos}\omega }_{1}tdt=\frac{1}{384}A\left(144a{A}^2+\left(192-24A^2+A^4\right)b-192\left(-1+{{\mathrm{\omega }}_2}^2\right)\right).$

(30)

He’s frequency formulation method provides

${\mathrm{\omega }}^2=\frac{{{\mathrm{\omega }}_1}^2\hspace{0.17em}\stackrel{\sim }{R_2}-{{\mathrm{\omega }}_2}^2\hspace{0.17em}\stackrel{\sim }{R_1}}{\stackrel{\sim }{R_2}-\stackrel{\sim }{R_1}}=1+\frac{3a{A}^2}{4}+b-\frac{A^{2}b}{8}+\frac{A^{4}b}{192}.$

(31)

Hence frequency of the Oscillator equation (15) according to He’s frequency formulation method is

$\omega =\sqrt{1+\frac{3a{A}^2}{4}+b-\frac{A^{2}b}{8}+\frac{A^{4}b}{192}}.$

(32)

Solution of equation (15) is

$x=A\cos \left(\sqrt{1+\frac{3a{A}^2}{4}+b-\frac{A^{2}b}{8}+\frac{A^{4}b}{192}}\right)t.$

(33)

Results and discussion

To demonstrate the efficiency of the present method solution Tables 1–3 are presented for $a=1,b=1$ and $A=2,\frac{\pi }{6},1$ , respectively. The columns of the tables represented exact solution (Runge–Kutta Fourth order), 2nd order approximate solution obtained by the MEBM, ³² 1st order approximate solution of the present method, 2nd order approximate solution obtained by the GRHB ³⁰ method and He’s frequency formulation solution along with the absolute percentage error sequentially. All the results of the table portrayed the supremacy of the present method over the other methods. Also, present results are fitted well with numerical results. Table 4 represents the nomenclature and we have considered here

$error\%=\left|\frac{Exact\hspace{0.17em}Value\hspace{0.17em}\left(\mathrm{R}\mathrm{u}\mathrm{n}\mathrm{g}\mathrm{e}-\mathrm{K}\mathrm{u}\mathrm{t}\mathrm{t}\mathrm{a}\hspace{0.17em}4\mathrm{t}\mathrm{h}\hspace{0.17em}\mathrm{o}\mathrm{r}\mathrm{d}\mathrm{e}\mathrm{r}\hspace{0.17em}\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\right)-Present\hspace{0.17em}value}{Exact\hspace{0.17em}Value\hspace{0.17em}\left(\mathrm{R}\mathrm{u}\mathrm{n}\mathrm{g}\mathrm{e}-\mathrm{K}\mathrm{u}\mathrm{t}\mathrm{t}\mathrm{a}\hspace{0.17em}4\mathrm{t}\mathrm{h}\hspace{0.17em}\mathrm{o}\mathrm{r}\mathrm{d}\mathrm{e}\mathrm{r}\hspace{0.17em}\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\right)}\right|\times 100.$

Table 1.

Comparison between the present method and other existing methods solution for $a=1,b=1$ and $A=2$ .

t	Exact x(t)	MEBM (2nd approx.) 32 error %	Present (1st approx.) error %	GRHB (2nd approx.) 30 error %	HFF (Equation (33)) error %
0.0	2.000000	2.000000 0.0000	2.000000 0.0000	2.000000 0.0000	2.000000 0.0000
0.5	0.9025057	.904695 0.24258	.9034881 0.108852	0.9015278 0.108354	0.959483 6.31323
1.0	−0.978461	−0.946665 3.249593	−0.9825991 0.422919	−0.981993 0.360975	−1.07939 10.3153
1.5	−1.9975295	−1.999243 0.085781	−1.9973918 0.006894	−1.9972898 0.012	−1.99514 0.119523
2.0	−0.8254517	−0.862431 4.479887	−0.8232094 0.271645	−0.819876 0.675473	−0.83491 1.14632
2.5	1.053172	0.988317 6.158063	1.06036579 0.683059	1.061086 0.751444	1.194056 13.3771
3.0	1.990132	1.996974 0.343796	1.9895809 0.027692	1.9891743 0.048122	1.98059 0.479481
3.5	0.747435	0.819899 9.695024	0.7419258 0.737081	0.737208 1.36828	0.706288 5.505107
4.0	−1.126481	−1.029620 8.598547	−1.1365982 0.898124	−1.138605 1.076272	−1.30292 15.6627
4.5	−1.977850	−1.993198 −.775994	−1.9766088 0.062755	−1.975698 0.108805	−1.95642 1.083671
5.0	−0.668582	−0.777122 16.23436	−0.6597869 1.315486	−0.653681 2.228747	−0.57423 14.11226

where error % denotes the absolute percentage error.

Table 2.

Comparison between the present method and other existing methods solution for $=1,$ $b=1$ , and $A=\frac{\pi }{6}$ .

t	Exact x(t)	MEBM (2nd approx.) 32 error %	Present (1st approx.) error %	GRHB (2nd approx.) 30 error %	HFF (Equation (33)) error %
0.0	0.523599	0.523599 0.0000	0.523599 0.0000	0.523599 0.0000	0.523599 0.0000
0.5	0.3854884	0.385491 0.000674	0.3854997 0.002931	0.38549851 0.002623	0.387775 0.59314
1.0	0.050207	0.050267 0.119505	0.0502051 0.003784	0.0502445 0.074691	0.050769 1.11957
1.5	−0.309795	−0.309673 0.039381	−0.3098009 0.001904	−0.3098001 0.001646	−0.31258 0.89766
2.0	−0.513458	−0.513428 0.005843	−0.513461 0.000584	−0.51346104 0.000592	−0.51375 0.05743
2.5	−0.4469344	−0.447019 0.018929	−0.446946 0.002595	−0.4469449 0.002349	−0.44839 0.32569
3.0	−0.148869	−0.149048 0.12024	−0.1488656 0.002284	−0.1488638 0.003493	−0.1504 1.02708
3.5	0.2228822	0.222615 0.119884	0.22288193 0.000121	0.22288195 0.000112	0.225623 1.22971
4.0	0.4835057	0.483386 0.024757	0.48351443 0.001806	0.48351456 0.001832	0.484588 0.22384
4.5	0.49153889	0.491650 0.022605	0.491546 0.001446	0.49154512 0.001267	0.492144 0.1231
5.0	0.24229490	0.242582 0.118492	0.24229483 0.0000289	0.24229215 0.001135	0.244371 0.85685

where error % denotes the absolute percentage error.

Table 3.

Comparison between the present method and other existing methods solution for $a=1,b=1$ and $A=1$ .

t	Exact x(t)	MEBM (2nd approx.) 32 error %	Present (1st approx.) error %	GRHB (2nd approx.) 30 error %	HFF (Equation (33)) error %
0.0	1.000000	1.000000 0.0000	1.000000 0.0000	1.000000 0.0000	1.00000 0.0000
0.5	0.6749066	0.674955 0.007171	0.6751213 0.03181	0.675057 0.02228	0.688849 2.06583
1.0	−0.047039	−0.045466 3.344034	−0.047043 0.0085	−0.047058 0.04039	−0.05097 8.36304
1.5	−0.743966	−0.742016 0.262109	−0.7442332 0.03592	−0.7442073 0.03243	−0.75908 2.0309
2.0	−0.9948065	−0.995120 0.031514	−0.9948181 0.00117	−0.9948117 0.00052	−0.9948 0.000352
2.5	−0.599965	−0.602381 0.40269	−0.6001034 0.02307	−0.5999978 0.00547	−0.61146 1.91666
3.0	0.140772	0.136096 3.321683	0.1407896 0.0125	0.140835 0.04475	0.15239 8.25306
3.5	0.8063175	0.802801 0.436119	0.8066064 0.03583	0.806614 0.03677	0.821412 1.87203
4.0	0.979312	0.980557 0.12713	0.979355 0.00439	0.9793295 0.00179	0.979268 0.004503
4.5	0.519968	0.525065 0.980253	0.5200191 0.00983	0.519871 0.01865	0.527724 1.49155
5.0	−0.233465	−0.225811 3.27844	−0.2335091 0.01889	−0.2335847 0.05127	−0.25222 8.03479

where error % denotes the absolute percentage error.

Table 4.

Nomenclatures.

Nomenclature	The meaning
$A$	Amplitude of the oscillator
$a$	Coefficient of the cubic nonlinear term
$b$	Coefficient of the harmonic restoring force
$x\left(t\right)$	Approximate solution
${\omega }_0$	Natural frequency
$\omega$	Frequency of the oscillator
$u,v,\cdots$	Unknown parameters for the MHBM
$m_1\left(\omega ,u,v,\cdots \right),m_3\left(\omega ,u,v,\cdots \right),\cdots$	Fourier series coefficients
HBM	Harmonic balance method
MEBM	Modified energy balance method
GRHB	Global residue harmonic balance method
MHBM	Modified harmonic balance method
HFF	He’s frequency formulation
$t$	Time
${\omega }_e$	Frequency obtained by the Runge-Kutta 4th order method for the Duffing oscillator
${\omega }_{Present}$	Frequency obtained by the present method for the Duffing oscillator
${\omega }_{HFF}$	Frequency obtained by He’s frequency formulation method for the Duffing oscillator

Absolute percentage error of the present method solution and the GRHB ³⁰ method solution are plotted in Figures 2–5 for $a=1,b=1$ and $A=\pi ,\frac{2\pi }{3},\frac{\pi }{2},\frac{\pi }{3}$ , respectively. Solution of the present method, numerical solution obtained by the Runge–Kutta 4th order method and the GRHB ³⁰ method solution for $a=1,b=1$ and $A=\pi$ are shown in Figures 6 and 7 for different time values.

Figure 2.

Comparison between the absolute percentage error obtained by the present method and the GRHB ³⁰ method solution for $a=1,$ $b=1$ ,and $A=\pi$ .

Figure 3.

Comparison between the absolute percentage error obtained by the present method and the GRHB ³⁰ method solution for $=\mathbf{1},$ $\mathit{b}=\mathbf{1}$ , and $\mathit{A}=\frac{\mathbf{2}\mathit{\pi }}{\mathbf{3}}$ .

Figure 4.

Comparison between the absolute percentage error obtained by the present method and the GRHB ³⁰ method for $\mathrm{a}=1,$ $\mathrm{b}=1$ and $\mathrm{A}=\frac{\mathrm{\pi }}{2}$ .

Figure 5.

Comparison between the absolute percentage error obtained by the present method and the GRHB ³⁰ method for $\mathrm{a}=1,$ $\mathrm{b}=1$ , and $\mathrm{A}=\frac{\mathrm{\pi }}{3}$ .

Figure 6.

Comparison among the solutions obtained by the present method, numerical method and the GRHB ³⁰ method, when $=1,$ $b=1$ , $A=\pi$ for $30\le t\le 32.5$ .

Figure 7.

Comparison among the solutions obtained by the present method, numerical method and the GRHB ³⁰ method when $=1,$ $b=1$ ,and $A=\pi$ for $5.4\le t\le 7$ .

Also, solution of the present method, numerical solution obtained by Runge–Kutta fourth order method and the GRHB ³⁰ method solution for $a=1,b=1$ and $A=\frac{\pi }{2},\frac{2\pi }{3},\frac{\pi }{3}$ are shown in Figures 8–10, respectively. From all the figures, it can be seen that only the first approximate solution of the present method provides better result than the second approximation solution of the GRHB ³⁰ method.

Figure 8.

Comparison among the solutions obtained by the present method, numerical method and the GRHB ³⁰ method, when $=1,$ $b=1$ , $A=\frac{2\pi }{3}$ for $30.2\le t\le 31.8$ .

Figure 9.

Comparison between solution obtained by the present method, numerical method and the GRHB ³⁰ method, when $=1,$ $b=1$ , and $A=\frac{\pi }{2}$ for $10\le t\le 15$ .

Figure 10.

Comparison among the solutions obtained by the present method, numerical method and the GRHB ³⁰ method when, $=1,$ $b=1$ , and $A=\frac{\pi }{3}$ for $4\le t\le 8$ .

In particular, when we put $a=1$ and $b=0$ in equation (15) then it reduces to the well-known Duffing oscillator^31,41,45

$\ddot{x}+x+x^3=0,$

(34)

The first approximation solution of equation (27) is

$x\left(t\right)=A\left(\left(1-u\right)\cos \phi +u\hspace{0.17em}\cos 3\phi \right),$

(35)

where $u$ and $\phi =\omega \hspace{0.17em}t$ is given by Equations (24) and (25), respectively, when $a=1$ and $b=0$ .

Figure 11 and Table 5 shows that the present method provides satisfactorily results for the Duffing oscillator for higher values of the amplitude when compared with the numerical results obtained by Runge–Kutta fourth order method and He’s frequency formulation method.

Figure 11.

Comparison among the solution obtained by the present method, He’s frequency formulation method and numerical method for Duffing oscillator, when $A=1000$ for $0\le t\le 0.0075$ .

Table 5.

Comparison of the frequency obtained by the present method, the numerical method and He’s frequency formulation method for different values of the amplitude $A$ for the Duffing oscillator.

A	${\omega }_e$	${\omega }_{Present}$ error (%)	${\omega }_{\mathrm{H}\mathrm{F}\mathrm{F}}$ (Equation (32)) error %
0.01	1.00004	1.00004 0.00000	1.00004 0.00000
0.1	1.00374	1.00374 0.00000	1.00374 0.00000
1.0	1.31778	1.31768 0.00714	1.32288 0.38702
5	4.35746	4.35495 0.05763	4.4441 1.98831
10	8.53359	8.53063 0.03462	8.7178 2.15865
100	84.72748	84.74443 0.02000	86.6083 2.21985
1000	847.21370	847.44964 0.02785	866.026 2.22049

where error % denotes the absolute percentage error.

Conclusion

A strongly nonlinear oscillator with cubic nonlinearity and with harmonic restoring force is solved by using a simple modified form of the harmonic balance method. The determination process of the present solution scheme is simple and straightforward compared to the modified energy balance method (MEBM), ³² the harmonic balance method, ¹⁵ and the GRHB ³⁰ method. Error analysis of the present method, the MEBM, ³² and the GRHB ³⁰ method was investigated and found decent results for the present method. For the present problem, only the first approximate solution shows a remarkable agreement when compared with the exact solution. Also, we have compared the present solution scheme with He’s frequency formulation solution and found superiority of the present method. Based on above determinations, the present method can be used to solve other strongly nonlinear oscillators effectively. The simple computational form and accuracy of the present method is obvious but the present method cannot solve the Van der Pol equation, the Mulholland equation etc.