Solving nonlinear differential equations of Vanderpol, Rayleigh and Duffing by AGM

Abstract

In the present paper, three complicated nonlinear differential equations in the field of vibration, which are Vanderpol, Rayleigh and Duffing equations, have been analyzed and solved completely by Algebraic Method (AGM). Investigating this kind of equations is a very hard task to do and the obtained solution is not accurate and reliable. This issue will be emerged after comparing the achieved solutions by numerical method (Runge-Kutte 4th). Based on the comparisons which have been made between the gained solutions by AGM and numerical method, it is possible to indicate that AGM can be successfully applied for various differential equations particularly for difficult ones. The results reveal that this method is not only very effective and simple, but also reliable, and can be applied for other complicated nonlinear problems.

Key words： Algebraic Method (AGM); Angular Frequency; Vanderpol; Rayleigh; Duffing

Introduction

Along with the rapid progress of nonlinear sciences, an intensifying interest among scientists and researchers has been emerged in the field of analytical asymptotic techniques particularly for nonlinear problems in the field of vibrations because this issue is very applicable in dynamics of structures (Mechanical, Earthquake, Civil Engineering take for example by Chopra [1]) and also in Electronic circuits in electrical engineering [2]. Although finding the solutions of linear equations by means of computer is very convenient, it is still very difficult and a time-consuming procedure to solve nonlinear problems either numerically or theoretically. Perhaps this is related to the fact that the various discredited methods or numerical simulations apply iteration techniques to find their numerical solutions of nonlinear problems and nearly all iterative methods are sensitive to initial solutions, so it is very difficult to obtain converged results in cases of strong nonlinearity. In addition, the most important information such as the natural circular frequency of a nonlinear oscillation depends on the initial conditions (i.e., amplitude of oscillation) will be lost during the procedure of numerical simulation. Perturbation methods [3,4] provide the most versatile tools available in nonlinear analysis of engineering problems and they are constantly being developed and applied to ever more complex problems. But like other nonlinear asymptotic techniques, perturbation methods have their own limitations take for example almost all perturbation methods are based on such an assumption that a small parameter must exist in an equation. This so-called small parameter assumption greatly restricts applications of perturbation techniques, as is well known, an overwhelming majority of nonlinear problems, especially those having strong nonlinearity, have no small parameters at all and so on.

Based on the above explanations, we should introduce some new developed methods for solving complicated nonlinear problems in different fields of study particularly in vibrations, where traditional techniques have not been successful up to now.

Furthermore, some techniques like perturbation methods are not practical for strongly nonlinear equations. As a result, due to conquer these weak-points, in recent years, much attention has been devoted to the newly developed manners to gain an approximate solution of nonlinear equations, such as Energy Balance Method [5,6], Homotopy Analysis Method [7,8], He’s Amplitude Frequency Formulation Method (HAFF) [9–14], Parameter-Expansion Method [15], Exp-function Method [10–19], Differential Transformation Method (DTM) [12,13], Homotopy Perturbation Method [11,16], Variational Iteration Method by J. H. He [14–17] and Adomian Decomposition Method [20,21] But the afore-mentioned methods do not have this ability to gain the solution of the presented problem in high precision. Therefore, these complicated nonlinear equations such as the presented problems in this paper should be solved by utilizing other approaches like AGM.

Analytical method

In general, vibrational equations and their initial conditions are defined for different systems as follows:

(1)f\left(\ddot{u},\dot{u},u,F_{\mathrm{0}}\sin ({\omega }_{\mathrm{0}}t)\right)=\mathrm{0},

(2)\left\{u(\mathrm{0})=A,\dot{u}(\mathrm{0})=\mathrm{0}\right\}.

Choosing the answer of the governing equation for solving differential equations by AGM

In AGM, a total answer with constant coefficients is required in order to solve differential equations in various fields of study, such as vibrations, structures, fluids and heat transfer. In vibrational systems with respect to the kind of vibration, it is necessary to choose the mentioned answer in AGM. To clarify here, we divide vibrational systems into two general forms:

1) Vibrational systems without any external force

Differential equations governing on this kind of vibrational systems are introduced in the following form:

(3)f\left(\ddot{u},\dot{u},u\right)=\mathrm{0.}

Now, the answer of this kind of vibrational system is chosen as

(4)u(t)=e^{-bt}\{A\cos (\omega t)+B\sin (\omega t)\}.

According to trigonometric relationships, Eq. (4) is rewritten as follows:

(5)u(t)={\mathrm{e}}^{-bt}\{a\cos (\omega t+\phi )\}.

It is notable that in the above equation a=\sqrt{A^{\mathrm{2}}+B^{\mathrm{2}}} and \phi =\arctan \left({\displaystyle \frac{B}{A}}\right).

Sometimes for increasing the precision of the considered answer of Eq. (3), we are able to add another term in the form of cosine by inspiration of Fourier cosine series expansion as follows:

(6)u(t)=e^{-bt}\{a\cos (\omega t+\phi )+d\cos (\mathrm{2}\omega t+\phi )\}.

In the above equation, we are able to omit the term (e^{-bt}) to facilitate the computational operations in AGM if the system is considered without any damping components.

Generally speaking in AGM, Eq. (5) or Eq. (6) is assumed as the answer of the vibrational differential Eq. (3) that its constant coefficients which are a, b, c \omega (angular frequency) and \phi (initial vibrational phase) can easily be obtained by applying the given initial conditions in Eq. (2). And also the above procedure will completely be explained through the presented example in the foregoing part of the paper.

It is noteworthy that if there is no damping component in the vibrational system, the constant coefficient b in Eqs. (5) and (6) will automatically be computed zero in AGM solution procedure.

On the contrary, the parameter b in Eqs. (5) and (6) for the other kind of vibrational system with damping component is obtained as a nonzero parameter in AGM.

2) Vibrational systems with external force

In this step, it is assumed that the external forces exerting on the vibrational systems are defined as

(7)F(t)=F_{\mathrm{0}}\sin ({\omega }_{\mathrm{0}}t).

As a result, the differential equation governing on the vibrational system is expressed like Eq. (1) as follows:

(8)f\left(\ddot{\mathit{u}},\dot{u},u,F_{\mathrm{0}}\sin ({\omega }_{\mathrm{0}}t)\right)=\mathrm{0.}

The answer of the above equation is introduced as the sum of the particular solution (u_p) and the harmonic solution (u_h) as follows:

(9)\begin{array}{l}{u}_h(t)=e^{-bt}\{A\cos (\omega t)+B\sin (\omega t)\},\hfill \\ u_p(t)=M\cos ({\omega }_{\mathrm{0}}t)+N\sin ({\omega }_{\mathrm{0}}t).\hfill \end{array}

Then

(10)u(t)=u_p+u_h.

By utilizing trigonometric relationships and substituting the yielded equations into Eq. (10), the desired answer will be obtained in the form of

(11)u(t)=e^{-bt}\{a\cos (\omega t+\phi )\}+d\cos ({\omega }_{\mathrm{0}}t+\varphi ).

To increase the precision of the achieved equation, we are able to add another term in the form of cosine by inspiration of Fourier cosine series expansion as follows:

(12)u(t)=e^{-bt}\{a\cos (\omega t+\phi )+c\cos (\mathrm{2}\omega t+\phi )\}+d\cos ({\omega }_{\mathrm{0}}t+\varphi ).

Finally, the exact solution of the all vibrational differential equations can be obtained in accordance with the following equation:

(13)u(t)=e^{-at}\left\{{\displaystyle \sum _{k=\mathrm{1}}^{\infty }{b}_k\cos (k\omega t+{\phi }_k)}\right\}+d\cos ({\omega }_{\mathrm{0}}t+\varphi ).

To deeply understand the above procedure, reading the following lines is recommended.

Since the constant coefficient (b) in vibrational systems without damping components is always obtained zero, we can add the term (b.t) Instead of (e^{-bt}) in Eq. (12) to decrease computational operations in the following form:

(14)u(t)=(b.t)+a\cos (\omega t+\phi )+d\cos (\mathrm{2}\omega t+\phi ).

Based on the above explanations, by applying initial conditions on a system without damping component, the value of parameter (b) is always zero for Eqs. (12) and (13). Therefore, the role of parameter (b) in both of Eqs.(12) and (13) which each of them can be considered as the answer of the vibrational problems is individually considered as a catalyst for increasing the precision of the assumed answer. However according to Eqs. (15) and (19) after applying initial conditions on the vibrational system in both states (with external force and without external force) by AGM, the value of parameter (b) is computed zero because the mentioned system has a free vibration without any damping component.

Again, we mention that in order to decrease computational operations for systems without damping components and since we know that (b) in the term (e^{-bt}) is zero so (e^{-bt}) can be omitted from Eq. (11). Consequently, Eq. (11) which has been considered as the answer of the systems without any damping component can be rewritten as follows:

(15)u(t)=a\cos (\omega t+\phi )+d\cos ({\omega }_{\mathrm{0}}t+\varphi ).

The constant coefficients of Eq. (11) or Eq. (12), which are \{a,b,\omega ,\phi ,c,d,\varphi \}, will easily be computed in AGM by applying the initial conditions of Eq. (2).

Application of initial conditions to compute constant coefficients and angular frequency by AGM

In AGM, the application of initial conditions of Eq. (2) is done in the two following forms:

a) Applying the initial conditions on the answer of differential equation

In regard to the kind of vibrational system (with external force and without external force) which was completely discussed in the previous part of this case study, a function is chosen as the answer of the differential equation from Eq. (5) or Eq. (6) for the systems without external forces and from Eq. (11) or Eq. (12) for the defined systems with external forces and then the initial conditions are applied on the selected function as follows:

(16)u(t)=u(IC).

It is notable that IC is the abbreviation of introduced initial conditions of Eq. (2).

b) Appling the initial conditions on the main differential equation and its derivatives

After choosing a function as the answer of differential equation according to the kind of vibrational system, this is the best time to substitute the mentioned answer into the main differential equation instead of its dependent variable (u).

Assume the general equation of the vibration such as Eq. (1) with time-independent parameter (t) and dependent function (u) as

(17)f\left(\ddot{\mathit{u}},\dot{\mathit{u}},u,F_{\mathrm{0}}\sin ({\omega }_{\mathrm{0}}t)\right)=\mathrm{0.}

Therefore, on the basis of the kind of vibrational system, a function as the answer of the differential equation, such as Eq. (5) or Eq. (6) and Eq. (11) or Eq. (12), are considered as follows:

(18)u=g(t).

In this step, the afore-mentioned equation is substituted into Eq. (16) instead of (u) in the following form:

(19)f(t)=f(g^{″}(t),g^{\prime }(t),g(t),F_{\mathrm{0}}\sin ({\omega }_{\mathrm{0}}t)).

Eventually, the application of initial conditions on Eq. (18) and its derivatives is expressed as

(20)\begin{array}{l}f(IC)=f(g^{″}(IC),g^{\prime }(IC),g(IC),\mathrm{...}),\hfill \\ f^{\prime }(IC)=f^{\prime }(g^{″}(IC),g^{\prime }(IC),g(IC),\mathrm{...}),\hfill \\ f^{″}(IC)=f^{″}(g^{″}(IC),g^{\prime }(IC),g(IC),\mathrm{...}),\hfill \\ \mathrm{...}\hfill \end{array}

To end up, it is better to say that in AGM after applying the initial conditions on Eq. (17), Eq. (18) and also on its derivatives from Eq. (19) according to the order of differential equation and utilizing the two given initial conditions of Eq. (2), a set of algebraic equations which is consisted of n equations with n unknowns is created. Therefore, the constant coefficients (a, b, c, d, angular frequency \omega and initial phase) are easily achieved, where this procedure will thoroughly be explained in the form of an example in the foregoing part of this paper.

It is noteworthy that in Eq. (19), we are able to use the derivatives of f(t) with higher orders until the number of yielded equations is equal to the number of the mentioned constant coefficients of the assumed answer.

Example 1

Consider the following nonlinear differential equation (Vanderpol, equation) in the form of

(21)f(t):{\displaystyle \frac{d^{\mathrm{2}}x}{d{t}^{\mathrm{2}}}}+ϵ(\mathrm{1}-x^{\mathrm{2}}=).{\displaystyle \frac{dx}{dt}}+{\beta }^{\mathrm{2}}x=\mathrm{0}

Then, the initial conditions are expressed as

(22)\mathit{x}(\mathrm{0})=A,\dot{x}(\mathrm{0})=\mathrm{0.}

It is notable that ϵ and β in the above equation are defined as constant values.

Solving the nonlinear differential equation by AGM

On the basis of the given explanations in the previous section (the analytical method), the answer of Eq. (21) is considered by AGM as a polynomials of Fourier series with constant coefficients as follows:

(23)x(t)=e^{-at}[b\cos (\omega t+\phi )].

It is notable that in AGM, the constant coefficients of Eq. (23), which are ω angular frequency), a, b, and \phi (initial vibrational phase), can easily be computed by applying initial or boundary conditions.

Applying initial or boundary conditions in AGM

Based on the given explanations in the previous section of this paper, the constant coefficients \omega ,a,b\mathrm{and}\phi from Eq. (23) are just achieved with respect to the given initial conditions and these initial conditions are applied in two manners in AGM.

a) With regard to Eq. (22), the initial conditions are applied on Eq. (23) as follows:

(24)\mathit{x}(\mathrm{0})=A,\mathrm{so}b\cos \phi =A.

Afterwards

(25)\dot{\mathit{x}}(\mathrm{0})=\mathrm{0},\mathrm{therefore}b(a\cos \phi +\omega \sin \phi )=\mathrm{0.}

b) In accordance with Eq. (20), the application of initial conditions on the main differential equation which in this case is Eq. (21) and also on its derivative is done after substituting Eq. (23) which has been considered as the answer of the main differential equation into Eq. (21) as

(26)$$\begin{gathered} f(\mathit{x}(\mathrm{0})): \\ (a^{\mathrm{2}}-{\omega }^{\mathrm{2}}+{\beta }^{\mathrm{2}}-\mathit{a}ϵ)\cos \phi +\omega (\mathrm{2}a-ϵ)\sin \phi +b^{\mathrm{2}}ϵ(\mathit{a}cos\phi +\omega sin\phi )co{s}^{\mathrm{2}}\phi =\mathrm{0.} \end{gathered}$$

Then for the first derivative of the achieved equation, we will have

(27)$$\begin{gathered} f^{\prime }(\mathit{x}(\mathrm{0})): \\ (a^{\mathrm{3}}-\mathrm{3}a{\omega }^{\mathrm{2}}-ϵa^{\mathrm{2}}+{\beta }^{\mathrm{2}}a+ϵ{\omega }^{\mathrm{2}})\cos \phi +\omega (\mathrm{3}{\mathit{a}}^{\mathit{2}}-{\omega }^{\mathrm{2}}-\mathrm{2}aϵ+{\beta }^{\mathrm{2}})\sin \phi \\ +b^{\mathrm{2}}\omega ϵ(\omega \sin \phi +\mathrm{3}a\cos \phi )\sin (\mathrm{2}\phi )+b^{\mathrm{2}}ϵ(\mathrm{3}{a}^{\mathrm{2}}-{\omega }^{\mathrm{2}}){\cos }^{\mathrm{3}}\phi =\mathrm{0.} \end{gathered}$$

By solving a set of algebraic equations which is consisted of four equations with four unknowns from Eqs. (24)-(27), the constant coefficients of Eq. (23) can easily be yielded.

For simplicity, the following new variables are considered:

(28)\psi =\mathrm{2}{ϵ}^{\mathrm{2}}{A}^{\mathrm{2}}+\mathrm{4}{\beta }^{\mathrm{2}}-{ϵ}^{\mathrm{2}}-{ϵ}^{\mathrm{2}}{A}^{\mathrm{4}}.

Then, the constant coefficients of Eq. (23) will be computed as

(29)a={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}ϵ(\mathrm{1}-A^{\mathrm{2}}),b={\displaystyle \frac{\mathrm{2}\beta A}{\sqrt{\psi }}},\phi =t{g}^{-\mathrm{1}}\left({\displaystyle \frac{\sqrt{\psi }}{\beta }}\right).

And then angular frequency is easily computed as

(30)\omega ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sqrt{\mathrm{2}{ϵ}^{\mathrm{2}}{A}^{\mathrm{2}}+\mathrm{4}{\beta }^{\mathrm{2}}-{ϵ}^{\mathrm{2}}-{ϵ}^{\mathrm{2}}{A}^{\mathrm{4}}}.

After substituting the obtained values from Eqs. (32) and (33) into Eq. (26), the answer of the presented problem is achieved in the following form:

(31)\mathit{x}(t)={\displaystyle \frac{\mathrm{2}\beta A}{\sqrt{\psi }}}{e}^{-\frac{\mathrm{1}}{\mathrm{2}}ϵ(\mathrm{1}-A^{\mathrm{2}})t}\cos \left\{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sqrt{\psi }t+t{g}^{-\mathrm{1}}\left({\displaystyle \frac{\sqrt{\psi }}{\beta }}\right)\right\}.

By selecting the following physical values:

(32)A=\mathrm{0.2},ϵ=\mathrm{0.15},\beta =\mathrm{1.2.}

So the answer and the related angular frequency of the presented problem, Eq. (21), are gained in the following form:

(33)\begin{array}{l}\omega =\mathrm{1.19784}\left({\displaystyle \frac{m}{s}}\right),\hfill \\ \mathit{x}(t)=\mathrm{0.20036}{e}^{-\mathrm{0.072}t}\cos (\mathrm{1.19784}t-\mathrm{0.060036}).\hfill \end{array}

According to the afore-mentioned equation, it is possible to draw the obtained solution in Figs. 1 and 2, and the related phase plane in Fig. 3 by AGM.

Fig.1 Chart of the obtained solution by AGM

Full size|PPT slide

Fig.2 Chart of the first derivative for the obtained solution by AGM

Full size|PPT slide

Fig.3 Resulted phase plane for example1 by AGM.

Full size|PPT slide

Numerical solution

In accordance with the given physical values from Eq. (32) and the introduced domain t\in \{\mathrm{0},\mathrm{50}\} which is defined in terms of second, the numerical solution (Runge-Kutte 4th) of the mentioned problem is presented in Table 1.

Tab.1 Obtained numerical solution of Eq. (31) based on the given physical values

t/s	0	10	20	30	40	50
u (t) Num. Rk 45	0.2	0.075603	0.0145714	-0.00544708	-0.0075291	-0.0047593
\dot{u}(t) Num. Rk 45	0.0	0.0632881	0.00493073	0.0249801	0.00849681	0.00109830

Comparing the obtained solutions by AGM and numerical method

Comparisons between the obtained solution by AGM and numerical method in Fig. 4 and the related derivative in Fig. 5 have been presented graphically and also the same procedure has been done for the relevant phase plane in Fig. 6.

Fig.4 Comparing the obtained solutions by AGM and numerical method

Full size|PPT slide

Fig.5 Comparing the first derivative of the obtained solutions by AGM and numerical method

Full size|PPT slide

Fig.6 Comparing the related phase planes of the achieved solutions by AGM and numerical method

Full size|PPT slide

In regard to the above charts, it is clear that AGM is a very applicable and reliable method for solving highly nonlinear vibrational differential equations with high precision like the presented problem in example 1.

Difference of the obtained solutions by AGM and numerical method

In Figs. 7 and 8, the differences between the obtained solution which is consisted of u(t)\mathrm{and}{u}^{\prime }(t) by AGM and numerical method have been presented.

Fig.7 Difference the obtained solutions by AGM and numerical method

Full size|PPT slide

Fig.8 Difference the first derivative of the obtained solutions by AGM and numerical method

Full size|PPT slide

Computation of other damping vibrational parameters

The maximum amplitude of vibrational displacement, velocity and acceleration equations can be computed as follows:

On the basis of the given explanations in Section 2.1, the answer of the introduced differential equation is considered as

(34)x(t)=b{e}^{-at}\{\cos (\omega t+\phi )\}.

In order to obtain the curve of locus for the maximum points of vibrational amplitudes we should consider that \cos (\omega t+\phi )=\pm \mathrm{1}, therefore we will have:

(35)x_{\max }(t)=b{\mathrm{e}}^{-at}.

In the above equation, the constant coefficients a and b can be obtained on the basis of the given explanations in Section 2.2. By substituting t_k={\displaystyle \frac{k\pi }{\omega }} in Eq. (35) in which k\in \{\mathrm{1},\mathrm{2},\mathrm{3},\cdots \}, it is possible to gain the maximum amplitudes such as x_{\mathrm{1}\max }, x_{\mathrm{2}\max } and x_{\mathrm{3}\max } According to the given physical values, the maximum amplitudes can be obtained as follows:

(36)x_{k,\max }=b{e}^{-a\left(\frac{k\pi }{\omega }\right)}.

Therefore, we will have

(37)k=\mathrm{1}\to t_{\mathrm{1}}={\displaystyle \frac{\pi }{\omega }},\mathrm{then}{x}_{\mathrm{1}\max }=b{e}^{-a\left(\frac{\pi }{\omega }\right)}=\mathrm{0.16588},

(38)k=\mathrm{2}\to t_{\mathrm{2}}={\displaystyle \frac{\mathrm{2}\pi }{\omega }},\mathrm{as}aresult,x_{\mathrm{2}\max }=b{e}^{-a\left(\frac{\mathrm{2}\pi }{\omega }\right)}=\mathrm{0.137338.}

The solution of the differential equation and the chart of locus for maximum displacements are shown simultaneously in Fig. 9.

Fig.9 Charts of the obtained solution and the locus of maximum displacement points

Full size|PPT slide

Furthermore, the maximum velocity {\dot{x}}_{\max } can be acquired after taking the first derivative of Eq. (34) as follows:

(39)\dot{\mathit{u}}(t)=b\{-a{e}^{-at}\cos (\omega t+\phi )-e^{-at}\omega \sin (\omega t+\phi )\}.

As regards trigonometric relations, the maximum velocity \left({\dot{x}}_{\max }(t)\right) from Eq. (39) can be obtained as follows:

(40){\dot{u}}_{\max }(t)=b\sqrt{a^{\mathrm{2}}+{\omega }^{\mathrm{2}}}{e}^{-at}.

In Eq. (40), t is defined as t_k={\displaystyle \frac{(\mathrm{2}k-\mathrm{1})\pi }{\mathrm{2}\omega }}, in which k\in \{\mathrm{1},\mathrm{2},\mathrm{3},\cdots \}. As a result, the equation of maximum points of vibrational velocity is expressed as follows:

(41)V={\dot{x}}_{k\max }=b\sqrt{a^{\mathrm{2}}+{\omega }^{\mathrm{2}}}{\mathrm{e}}^{-a\frac{(\mathrm{2}k-\mathrm{1})\mathrm{\pi }}{\mathrm{2}\omega }}.

It is necessary to mention that a, b and \omega have been obtained with respect to the given physical values. Consequently, the maximum vibrational velocities can be computed as follows:

(42)k=\mathrm{1}\to {\dot{x}}_{\mathrm{1}\max }=\mathrm{0.21877}\left({\displaystyle \frac{m}{\mathrm{s}}}\right),\mathrm{a}\mathrm{n}\mathrm{d}k=\mathrm{2}\to {\dot{x}}_{\mathrm{2}\max }=\mathrm{0.1811252}\left({\displaystyle \frac{m}{\mathrm{s}}}\right).

Then, the chart of locus for vibrational velocity and the own vibrational velocity chart are illustrated in Fig. 10.

Fig.10 Achieved results for the vibrational velocity and the locus of the maximum vibrational velocity points

Full size|PPT slide

Moreover, the desired equation for vibrational acceleration will be obtained by computing the second derivative of the obtained solution or the first derivative of the velocity equation. And also by utilizing trigonometric equations, it is possible to compute the locus for maximum acceleration points in t_k={\displaystyle \frac{k\pi }{\omega }} as follows:The locus of maximum acceleration points

(43)\mathrm{T}\mathrm{h}\mathrm{e}\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{u}\mathrm{s}\mathrm{o}\mathrm{f}\max \mathrm{i}\mathrm{m}\mathrm{u}\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{p}\mathrm{o}\mathrm{int}\mathrm{s}\to \ddot{u}(t)=b(a^{\mathrm{2}}+{\omega }^{\mathrm{2}})e^{-at}.

The maximum acceleration points

(44)\mathrm{T}\mathrm{h}\mathrm{e}\max \mathrm{i}\mathrm{m}\mathrm{u}\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{p}\mathrm{o}\mathrm{int}\mathrm{s}\to {\ddot{u}}_{k\max }=b(a^{\mathrm{2}}+{\omega }^{\mathrm{2}})e^{-a\left(\frac{k\pi }{\omega }\right)}.

In regard to the given physical values, the maximum vibrational accelerations in damping states can be obtained as follows:

(45)k=\mathrm{1}\to {\ddot{x}}_{\mathrm{1}\max }=\mathrm{0.238872}\left({\displaystyle \frac{m}{{\mathrm{s}}^{\mathrm{2}}}}\right),\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{f}\mathrm{o}\mathrm{r}k=\mathrm{2}\to {\ddot{x}}_{\mathrm{2}\max }=\mathrm{0.197767}\left({\displaystyle \frac{m}{{\mathrm{s}}^{\mathrm{2}}}}\right).

Therefore, the charts of locus for maximum acceleration points and the obtained vibrational acceleration equation are depicted in Fig.11.

Fig.11 Resulted charts of the vibrational acceleration equation and its related locus

Full size|PPT slide

The damping factor in the vibrational solution of b\{e^{-at}\cos (\omega \mathit{t}+\phi )\} is e^{-at}, since in the vibrational differential equation, the term e^{-\xi \omega t} is the main factor for damping, (\xi ) can be obtained as follows:

(46)e^{-at}\equiv e^{-\xi \omega t}\mapsto \xi ={\displaystyle \frac{a}{\omega }}.

Consequently, according to Eqs. (29) and (30) damping ratio is calculated as follows:

(47)\xi ={\displaystyle \frac{ϵ(\mathrm{1}-A^{\mathrm{2}})}{\sqrt{\mathrm{2}{ϵ}^{\mathrm{2}}{A}^{\mathrm{2}}+\mathrm{4}{\beta }^{\mathrm{2}}-{ϵ}^{\mathrm{2}}-{ϵ}^{\mathrm{2}}{A}^{\mathrm{4}}}}}.

In Eq. (46), the parameters a and \omega have been gained in accordance with the given physical values, and the initial vibrational phase can be obtained as follows:

(48)\phi =t{g}^{-\mathrm{1}}(\xi ).

Therefore in this case study, it is clear that \xi =\mathrm{0.0601083}.

Results and discussion

The chart of the angular frequency (\omega ) from Eq. (30) in terms of initial amplitude of vibration in accordance with the introduced variables from Eq. (32) is illustrated in Fig.12.

Fig.12 Chart of angular frequency in terms of initial vibrational amplitude

Full size|PPT slide

It is clear that the more amount of amplitude of vibration in the initial condition (IC), the more increasing of the angular frequency.

The chart of damping ratio (\xi ) in terms of initial amplitude of vibration (A) from Eq. (47) is depicted in Fig.13.

Fig.13 Variation of damping ratio in terms of the initial amplitude of vibration (A)

Full size|PPT slide

With regard to Fig. 13, it is revealed that by increasing the amount of initial amplitude of vibration (A), the value of damping ratio will be decreased. Since, the damping coefficient is computed by C=\mathrm{2}m\omega \xi, there is a direct relationship between initial amplitude of vibration and damping coefficient.

Example 2

Consider the following nonlinear differential equation (Rayleigh equation) in the form of

(49)f(t):\ddot{x}+\left(\alpha -\beta {\dot{x}}^{\mathrm{2}}\right)\dot{x}+{\lambda }^{\mathrm{2}}x=\mathrm{0.}

Then, the initial conditions are expressed as

(50)\mathit{x}(\mathrm{0})=A,\dot{x}(\mathrm{0})=\mathrm{0.}

It is notable that \alpha, \beta, and \lambda in the above equation are defined as constant values.

Solving the differential equation by AGM

To solve Eq. (49), a finite series with constant coefficients has been considered as the answer of follows:

(51)x(t)=e^{-at}\{b\cos (\omega t+\phi )\}.

Constant coefficients a, b, \omega, \phi in the afore-mentioned equation will be computed by applying initial conditions.

Applying initial conditions in AGM

Exactly like the procedure which has been explained in example 1, the initial conditions are applied in two ways in AGM. Therefore, we have

(52)\mathit{x}(\mathrm{0})=A,\mathrm{so}b\cos \phi =\mathrm{0.}

Afterwards

(53)\dot{\mathit{x}}(\mathrm{0})=\mathrm{0},\mathrm{therefore}b(a\cos \phi +\omega \sin \phi )=\mathrm{0.}

Then, applying the initial conditions on the main differential equation which is Eq. (49) and is named f(\mathit{t}) and on its derivatives is done after substituting Eq. (51) into differential Eq. (49) in the following form:

1) applying initial conditions on the yielded equation:

(54)$$\begin{gathered} f(\mathit{x}(\mathit{t}=\mathrm{0})): \\ ({\mathit{a}}^{\mathrm{2}}-{\omega }^{\mathrm{2}}+{\lambda }^{\mathrm{2}}-\alpha \mathit{a})cos\phi +\omega (\mathrm{2}a-\alpha )\sin \phi \\ +{\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}a{b}^{\mathrm{2}}\beta \omega (acos\phi +\omega sin\phi )sin(\mathrm{2}\phi )+b^{\mathrm{2}}\beta (a^{\mathrm{3}}{\cos }^{\mathrm{3}}\phi +{\omega }^{\mathrm{3}}{\sin }^{\mathrm{3}}\phi )=\mathrm{0.} \end{gathered}$$

2) applying initial conditions on the first derivative of the yielded equation:

(55)$$\begin{gathered} f^{\prime }(x(t=\mathrm{0})): \\ (a^{\mathrm{3}}-\mathrm{3}a{\omega }^{\mathrm{2}}-\alpha a^{\mathrm{2}}+\alpha {\omega }^{\mathrm{2}}+{\lambda }^{\mathrm{2}}a)\cos \phi +\omega (\mathrm{3}{a}^{\mathrm{2}}-{\omega }^{\mathrm{2}}-\mathrm{2}\alpha a+{\lambda }^{\mathrm{2}})\sin \phi \\ +{\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}{b}^{\mathrm{2}}\beta \omega (\mathrm{5}{a}^{\mathrm{2}}\omega \sin \phi -{\omega }^{\mathrm{2}}\sin \phi -\mathrm{2}a{\omega }^{\mathrm{2}}\cos \phi +\mathrm{4}{a}^{\mathrm{3}}\cos \phi )\sin (\mathrm{2}\phi ) \\ +\mathrm{3}{a}^{\mathrm{2}}b{}^{\mathrm{2}}\beta (a^{\mathrm{2}}-{\omega }^{\mathrm{2}}){\cos }^{\mathrm{3}}\phi +\mathrm{6}\beta b^{\mathrm{2}}a{\omega }^{\mathrm{3}}{\sin }^{\mathrm{3}}\phi =\mathrm{0.} \end{gathered}$$

In this step, by solving a set of algebraic equations which is consisted of four equations with four unknowns from Eqs. (52) to (55), the constant coefficients of Eq. (51) which are a, b, \omega, \phi will be obtained very easily as follows:

(56)a={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\alpha ,b={\displaystyle \frac{\mathrm{2}A\lambda }{\sqrt{\mathrm{4}{\lambda }^{\mathrm{2}}-{\alpha }^{\mathrm{2}}}}},\omega ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sqrt{\mathrm{4}{\lambda }^{\mathrm{2}}-{\alpha }^{\mathrm{2}}},\phi =-t{g}^{-\mathrm{1}}\left({\displaystyle \frac{\alpha }{\lambda }}\right).

After substituting the obtained values from Eq. (56) into Eq. (51), the answer of the presented problem is achieved in the following form:

(57)x(t)={\displaystyle \frac{\mathrm{2}\lambda A}{\sqrt{\mathrm{4}{\lambda }^{\mathrm{2}}-{\alpha }^{\mathrm{2}}}}}{e}^{-\frac{\mathrm{1}}{\mathrm{2}}\alpha t}\cos \left\{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sqrt{\mathrm{4}{\lambda }^{\mathrm{2}}-{\alpha }^{\mathrm{2}}}t-t{g}^{-\mathrm{1}}\left({\displaystyle \frac{\alpha }{\lambda }}\right)\right\}.

By selecting the following physical values:

(58)\begin{array}{llll}A=\mathrm{0.1},\hfill & \alpha =\mathrm{0.2},\hfill & \beta =\mathrm{0.1},\hfill & \lambda =\mathrm{1.2}\hfill \end{array}.

So the answer and the related angular frequency of the presented problem, Eq. (49), are gained in the following form:

(59)\begin{array}{l}\omega =\mathrm{1.1958},\hfill \\ x(t)=\mathrm{0.10035}{e}^{-\mathrm{0.1}t}\cos (\mathrm{1.1958}t-\mathrm{0.08343}).\hfill \end{array}

Consequently, the charts of the obtained solution by AGM in Fig. 14 and its first derivative in Fig. 15 are depicted.

Fig.14 Chart of the obtained solution by AGM

Full size|PPT slide

Fig.15 Chart of the first derivative for the obtained solution by AGM

Full size|PPT slide

Afterwards, the related phase plane is depicted in Fig. 16.

Fig.16 Resulted phase plane for example 2 by AGM

Full size|PPT slide

Numerical solution of the differential equation

The achieved results of numerical solution are presented in the specified domain t\in \{\mathrm{0},\mathrm{40}\} in Table2.

Tab.2 Results of numerical solution based on the given physical values in the specified domain

t	0	8	16	24	32	40
u(t) Num. Rk 45	0.1	-0.045104	0.01990540	-0.00859812	0.0036300	-0.00149353
u^{\prime }(t) Num.Rk 45	0.0	0.0076792	-0.00662959	0.00452472	-0.00264644	0.001440892

Comparing the achieved solutions by AGM and numerical method

In this step, Figs. 17 and 18 are depicted according to Table 2 from numerical method and the obtained solution by AGM in order to compare the achieved solutions together.

Fig.17 Comparison between the achieved solutions by AGM and numerical method

Full size|PPT slide

Fig.18 Comparing the first derivative of the obtained by AGM and numerical method

Full size|PPT slide

After that, the obtained phase plane by AGM has been compared with numerical solution in Fig. 19.

Fig.19 Comparing the related phase planes of the achieved solutions by AGM and numerical method

Full size|PPT slide

Difference of the obtained solutions by AGM and numerical method

The existed differences between the yielded solution by AGM and numerical method are illustrated in Figs. 20 and 21.

Fig.20 Difference of the obtained solutions by AGM and numerical method

Full size|PPT slide

Fig.21 Difference of the first derivative of the obtained solutions by AGM and numerical method

Full size|PPT slide

Example 3

Consider the following nonlinear differential equation (Duffing equation) in the form of

(60)f(t):{\displaystyle \frac{d^{2}u}{d{t}^{\mathrm{2}}}}+{\beta }^{\mathrm{2}}(u+{\mu }^{\mathrm{2}}{u}^{\mathrm{3}})=P\sin (\mathrm{\Omega }t).

Then, the initial conditions are expressed as

(61)u(\mathrm{0})=A,\dot{u}(\mathrm{0})=\mathrm{0.}

It is notable that \beta, \mu, \mathrm{\Omega } and P in the above equation are defined as constant values are expressed as

(62)\begin{array}{lllll}A=\mathrm{0.15},\hfill & \mathrm{\Omega }=\mathrm{2},\hfill & \beta =\mathrm{0.3},\hfill & \mu =\mathrm{0.2}\hfill & P=\mathrm{0.2}\hfill \end{array}.

Solving the differential equation by AGM

To solve Eq. (60), a finite series with constant coefficients has been considered as the answer of follows:

(63)u(t)=e^{-at}\{b\cos (\omega t+\phi )\}+dsin(\mathrm{\Omega }t+\varphi ).

Constant coefficients a, b, d, \omega, \phi and \varphi in the afore-mentioned equation will be computed by applying initial conditions.

Applying initial conditions in AGM

Exactly like the procedure which has been explained in examples 1 and 2, the initial conditions are applied in two ways in AGM. Therefore, we have

(64)u(\mathrm{0})=A,\mathrm{so}b\cos \phi +d\sin \varphi =\mathrm{0.15.}

Afterwards

(65)\dot{u}(\mathrm{0})=\mathrm{0},\mathrm{therefore}-b(a\cos \phi +\omega \sin \phi )+\mathrm{2}dcos\varphi =\mathrm{0.}

Then, applying the initial conditions on the main differential equation which is Eq. (60) and is named f(t) and on its derivatives is done after substituting Eq. (63) into differential Eq. (60) in the following form:

1) applying initial conditions on the yielded equation:

(66)$$\begin{gathered} f(u(t=\mathrm{0})): \\ \mathrm{so} \\ (a^{\mathrm{2}}-{\omega }^{\mathrm{2}}+\mathrm{0.09})cos\phi +\mathrm{2}ba\omega \sin \phi -\mathrm{3.91}d\sin \varphi +\mathrm{0.0036}{(b\cos \phi +d\sin \varphi )}^{\mathrm{3}}=\mathrm{0.} \end{gathered}$$

2) applying initial conditions on the first derivative of the yielded equation:

(67)$$\begin{gathered} f^{\prime }(u(t=\mathrm{0})): \\ ba(\mathrm{3}{\omega }^{\mathrm{2}}-a^{\mathrm{2}}-\mathrm{0.09})\cos \phi +b({\omega }^{\mathrm{3}}-\mathrm{3}{a}^{\mathrm{2}}\omega -\mathrm{0.09}\omega ))\sin \phi -\mathrm{7.82}d\cos \varphi \\ -\mathrm{0.0108}{(b\cos \phi +d\sin \varphi )}^{\mathrm{2}}(ab\cos \phi +b\omega \sin \phi -\mathrm{2}d\cos \varphi )=\mathrm{0.4.} \end{gathered}$$

3) applying initial conditions on the second derivative of the yielded equation:

(68)$$\begin{gathered} f^{″}(u(t=\mathrm{0})): \\ b(a^{\mathrm{4}}-\mathrm{6}{a}^{\mathrm{2}}{\omega }^{\mathrm{2}}+{\omega }^{\mathrm{4}}-\mathrm{0.09}{\omega }^{\mathrm{2}}+\mathrm{0.09}{a}^{\mathrm{2}})\cos \phi +ab(\mathrm{4}{a}^{\mathrm{2}}\omega -\mathrm{4}{\omega }^{\mathrm{3}}+\mathrm{0.18}\omega )\sin \phi \\ +\mathrm{15.64}d\sin \varphi +\mathrm{0.0216}(b\cos \phi +d\sin \varphi ){(ab\cos \phi +b\omega \sin \phi -\mathrm{2}d\cos \varphi )}^{\mathrm{2}} \\ +\mathrm{0.0108}{(bcos\phi +dsin\varphi )}^{\mathrm{2}}(a^{\mathrm{2}}b\cos \phi +\mathrm{2}ab\omega \sin \phi -b{\omega }^{\mathrm{2}}\cos \phi -\mathrm{4}d\sin \varphi )=\mathrm{0.} \end{gathered}$$

4) applying initial conditions on the third derivative of the yielded equation:

(69)$$\begin{gathered} f^{‴}(u(t=\mathrm{0})): \\ ab(\mathrm{10}{a}^{\mathrm{2}}{\omega }^{\mathrm{2}}-a^{\mathrm{4}}-\mathrm{5}{\omega }^{\mathrm{4}}-\mathrm{0.09}{a}^{\mathrm{2}}+\mathrm{0.27}{\omega }^{\mathrm{2}})\cos \phi +b(\mathrm{10}{a}^{\mathrm{2}}{\omega }^{\mathrm{2}}-\mathrm{5}\omega a^{\mathrm{4}}-{\omega }^{\mathrm{5}}-\mathrm{0.27}{a}^{\mathrm{2}}\omega +\mathrm{0.09}{\omega }^{\mathrm{3}})\sin \phi \\ +\mathrm{31.28}d\cos \varphi -\mathrm{0.0216}{(ab\cos \phi +b\omega \sin \phi -\mathrm{2}d\cos \varphi )}^{\mathrm{3}}+\mathrm{0.0648}(b\cos \phi +d\sin \varphi )(-ab\cos \phi -b\omega \sin \phi +\mathrm{2}d\cos \varphi ) \\ \{b(a^{\mathrm{2}}-{\omega }^{\mathrm{2}})\cos \phi +\mathrm{2}ab\omega \sin \phi -\mathrm{4}d\sin \varphi \}+\mathrm{0.0108}{(b\cos \phi +d\sin \varphi )}^{\mathrm{2}}\{ab(\mathrm{3}{\omega }^{\mathrm{2}}-a^{\mathrm{2}})\cos \phi +b\omega ({\omega }^{\mathrm{2}}-\mathrm{3}{a}^{\mathrm{2}})\sin \phi -\mathrm{8}d\cos \varphi \}=-\mathrm{1.6.} \end{gathered}$$

In this step, by solving a set of algebraic equations which is consisted of six equations with six unknowns from Eqs. (64) to (69), the constant coefficients of Eq. (63) which are a, b, d, \omega, \phi and \varphiwill be obtained very easily as follows:

(70)a=-\mathrm{0.000101788},b=\mathrm{0.372137},d=\mathrm{0.051154},\omega =\mathrm{0.30036},\phi =-\mathrm{1.15592},\varphi =\mathrm{0.}

After substituting the obtained values from Eq. (70) into Eq. (63), the answer of the presented problem is achieved in the following form:

(71)u(t)=\mathrm{0.372137}{\mathrm{e}}^{\mathrm{0.000101788}t}\cos (\mathrm{0.30036}t-\mathrm{1.15592})-\mathrm{0.051154}\sin (\mathrm{2}t).

Consequently, the charts of the obtained solution and its first derivative are depicted in Figs. 22 and 23 respectively.

Fig.22 Chart of the obtained solution by AGM

Full size|PPT slide

Fig.23 Chart of the first derivative for the obtained solution by AGM

Full size|PPT slide

Then, the related phase plane has been depicted in Fig. 24.

Fig.24 The resulted phase plane for the presented Duffing equation by AGM.

Full size|PPT slide

Numerical solution of the differential equation

The achieved results of numerical solution are presented in the specified domain t\in \{\mathrm{0},\mathrm{40}\} in Table 3.

Tab.3 Results of numerical solution based on the given physical values in the specified domain

t	0	8	16	24	32	40
u(t)Num. Rk 45	0.15	0.132085	-0.3522582	0.40210454	-0.26044564	0.004382949
u^{\prime }(t)Num. Rk 45	0.00	-0.0081056	-0.0304086	0.09030328	-0.13168945	0.12219430

Comparing the achieved solutions by AGM and numerical method

The following charts which are consisted of Figs. 25 and 26 are depicted according to Table 3 from numerical method and the obtained solution by AGM in order to compare them together.

Fig.25 Comparison between the achieved solutions by AGM and numerical method

Full size|PPT slide

Fig.26 Comparing the first derivative of the obtained by AGM and numerical method

Full size|PPT slide

In this step, the yielded phase plane by AGM has been compared with numerical solution in Fig. 27.

Fig.27 Comparing the related phase planes of the achieved solutions by AGM and numerical method

Full size|PPT slide

Difference of the obtained solutions by AGM and numerical method

The existed differences between the yielded solution by AGM and numerical method based on Table 3 are illustrated in Figs. 28 and 29.

Fig.28 Difference the obtained solutions by AGM and numerical method

Full size|PPT slide

Fig.29 Difference the first derivative of the obtained solutions by AGM and numerical method

Full size|PPT slide

Conclusions

In this paper, three complicated nonlinear vibrational differential equations which are Vanderpol , Rayleigh and Duffing equations have been introduced and analyzed completely by Algebraic Method (AGM) and also the obtained results have been compared with numerical method. Then, the vibrational velocity and vibrational acceleration have successfully been achieved. Afterwards, the related equations of locus for vibrational velocity and acceleration have been gained and depicted completely. Eventually, the equation of damping ratio in terms of initial amplitude of vibration and angular frequency has been obtained perfectly. The above process has been done in order to show the ability of AGM for solving a broad range of differential equations in different fields of study particularly in vibrations. Consequently, it is concluded that AGM is a reliable and precise approach for solving miscellaneous differential equations. Moreover, a summary of the AGM excellence and benefits is explained as: By solving a set of algebraic equations with constant coefficients, we are able to obtain the solution of nonlinear differential equation along with the related angular frequency simultaneously very easily which applying this procedure is possible even for students with intermediate mathematical knowledge. On the other hand, it is better to say that AGM is able to solve linear and nonlinear differential equations directly in most of the situations that means the final solution can be obtained without any dimensionless procedure. Therefore, AGM can be considered as a significant progress in nonlinear sciences.

Acknowledgements

The authors are grateful to the Ancient Persian mathematician, astronomer and geographer Muammadibn Musa Kharazmi, who was the one that his Compendious Book on Calculation by Completion and Balancing presented the first systematic solution of linear and quadratic equations. Furthermore, he was considered as the original inventor of algebra and is the one who Europeans derive the term ALGEBRA from his book and also the expression ALGORITHM has been taken from his name (the Latin form of his name). Consequently, we dedicate this way of solving linear and nonlinear differential equations to this scientist.

References

1 ChopraA K. Dynamic of Structures, Theory and Applications to Earthquake Engineering. New Jersey: Prentice-Hall, 1995

2 ZhouJ K. Differential Transformation and Its Applications for Electrical Circuits. Wuhan: Huazhong University Press, 1986, 363–364

3 ActonJ R, SquireP T. Solving Equations with Physical Understanding. Boston: Adam Hilger Ltd, 1985

4 BellmanR. Perturbation Techniques in Mathematics, Physics and Engineering. New York: Holt, Rinehart and Winston, Inc, 1964, 254–267

5 SfahaniM G, BarariA, OmidvarM, GanjiS S, DomairryG. Dynamic response of inextensible beams by improved energy balance method. Proceedings of the Institution of Mechanical Engineers, Part K: Journal of multi-body Dynamics, 2011, 225(1): 66–73

6 MomeniM, JamshidiN, BarariA, GanjiD D. Application of He’s energy balance method to duffing harmonic oscillators. International Journal of Computer Mathematics, 2011, 88(1): 135–144

7 GanjiZ Z, GanjiD D, BararniaH. Approximate general and explicit solutions of nonlinear BBMB equations exp-function method. Applied Mathematical Modelling, 2009, 33(4): 1836–1841

8 PirbodaghiT, AhmadianM, FesangharyM. On the homotopy analysis method for non-linear vibration of beams. Mechanics Research Communications, 2009, 36(2): 143–148

9 HeJ H. An improved amplitude-frequency formulation for nonlinear oscillators. International Journal of Nonlinear Sciences and Numerical Simulation, 2008, 9(2): 211–223

10 GanjiS S, GanjiD D, SfahaniM G, KarimpourS. Application of HAFF and HPM to the systems of strongly nonlinear oscillation. Current Applied Physics, 2010, 10(5): 1317–1325

11 HeJ H. The homotopy perturbation method for nonlinear oscillators with discontinuities. Applied Mathematics and Computation, 2004, 151: 287–292

12 HeJ H, WuX H. Exp-function method for nonlinear wave equations. Chaos. Solitos & Fractals, 2006, 30(3): 700–708

13 ChangJ R. The exp-function method and generalized solitary solutions. Computers & Mathematics with Applications, 2011, 61(8): 2081–2084

14 GanjiS S, GanjiD D, BabazadehH, SadoughiN. Application of amplitude-frequency formulation to nonlinear oscillation system of the motion of a rigid rod rocking back. Mathematical Methods in the Applied Sciences, 2010, 33(2): 157–166

15 GanjiD D. The application of He’s homotopy perturbation method to nonlinear equations arising in heat transfer. Physics Letters, 2006, 355(4-5): 337–341

16 HeJ H. Homotopy perturbation technique. Computer Methods in Applied Mechanics and Engineering, 1999, 178(3): 257–262

17 GhoshS, RoyA, RoyD. An adaptation of adomian decomposition for numeric-analytic integration of strongly nonlinear and chaotic oscillators. Computer Methods in Applied Mechanics and Engineering, 2007, 196(4–6): 1133–1153

18 HeJ H. Approximate analytical solution for seepage flow with fractional derivatives in porous media. Computer Methods in Applied Mechanics and Engineering, 1998, 167(1–2): 57–68

19 HeJ H. Approximate analytical solution of Blasiu’s equation. Communications in Nonlinear Science and Numerical Simulation, 1998, 3(4): 260–263

20 AbbasbandyS. A numerical solution of Blasius equation by Adomian’s decomposition method and comparison with homotopy perturbation method. Chaos, Solitons, and Fractals, 2007, 31(1): 257–260

21 AdomianG. Solving Frontier Problems of Physics: The Decomposition Method. Boston, MA: Kluwer Academic Publishers, 1994