The simpler,the better: Analytical methods for nonlinear oscillators and fractional oscillators

Abstract

In engineering, a fast estimation of the periodic property of a nonlinear oscillator is much needed. This paper reviews some simplest methods for nonlinear oscillators, including He’s frequency formulation, the max-min approach and the homotopy perturbation method. A mathematical insight into He’s frequency formulation is given, and the weighted average is introduced to further improve the estimated accuracy of the frequency. Fractional oscillators are also discussed.

Introduction

Nonlinear oscillations arise everywhere in our everyday life and engineering. As an exact solution might be too complex to be used for a practical application, many analytical methods have been used in open literature, for example, the variational iteration method,^1–7 the homotopy perturbation method,^8–20 He–Laplace method,^21–23 the variational approach^24–29 and the Hamiltonian approach.^30,31 The most important property of a nonlinear oscillator is the relationship between the frequency and its amplitude, the simplest method to estimate the frequency–amplitude relationship might be He’s frequency formulation^32–34 and the max-min approach,^35,36 which are still under development and many modifications were proposed to improve the accuracy.^37–45

Consider a nonlinear oscillator in the form $$u^{″}+f(u)=0\text{,\hspace{1em}}u(0)=A,u^{\prime }(0)=0$$ (1)

The simplest criterion for the existence of a periodic solution of equation (1) is ¹ $$\frac{f(u)}{u}>0$$ (2)

This criterion can be understood by a pendulum-like oscillator $$u^{″}=-\left\{\frac{f(u)}{u}\right\}u$$ (3)where $u^{″}$ is the acceleration, u is the displacement. When the acceleration is positive/negative, its displacement must be negative/positive, otherwise the displacement might tend to be either infinity or zero, and no periodic motion might occur. As an example, we consider the following equation $$u^{″}+\mathrm{sign}(u)=0\text{,\hspace{1em}}u(0)=A,u^{\prime }(0)=0$$ (4)where $\mathrm{sign}(u)$ is a sign function defined as $$\mathrm{sign}(u)=\{\begin{array}{l}1,u>0\hfill \\ 0,u=0\hfill \\ -1,u<0\hfill \end{array}$$ (5)It is obvious that we have $$\frac{\mathrm{sign}(u)}{u}>0$$ (6)

So equation (4) has a periodic solution.

In this paper, we will review some simplest methods for fast calculation of the frequency–amplitude relationship.

He’s frequency formulation

He’s frequency formulation was initially inspired by an ancient Chinese algorithm ³⁶ ; now this formulation has gained much attention and has been under development, with various modifications appearing in open literature.^37–45

We define a residual function R(t) as $$R(t)=u^{″}+f(u)$$ (7)

He’s frequency formulation is to choose two trial solutions $$u_1\text{(}t\text{)}=A\cos {\omega }_{1}t$$ (8) $$u_2\text{(}t\text{)}=A\cos {\omega }_{2}t$$ (9)where ${\omega }_1$ and ${\omega }_2$ are the arbitrarily chosen frequencies. The residuals are, respectively, as $$R_1(t)=-A{\omega }_1^2\cos {\omega }_{1}t+f(A\cos {\omega }_{1}t)$$ (10) $$R_2(t)=-A{\omega }_2^2\cos {\omega }_{2}t+f(A\cos {\omega }_{2}t)$$ (11)

The residuals are time-dependent, and we can use their weighted average $${\tilde{R}}_1={\displaystyle {\int }_0^{T/4}{R}_1(t)w(t)\text{d}t}$$ (12) $${\tilde{R}}_2={\displaystyle {\int }_0^{T/4}{R}_2(t)w(t)\text{d}t}$$ (13)where w(t) is a weighted function. He’s frequency formulation is $${\omega }^2=\frac{{\tilde{R}}_2{\omega }_1^2-{\tilde{R}}_1{\omega }_2^2}{{\tilde{R}}_2-{\tilde{R}}_1}$$ (14)

If the weighted function is chosen as $$w(t)=\cos \omega t$$ (15)we have the original He’s formulation. If we choose w(t)=1, we have Ren–Hu’s modification. ⁴¹ If we choose $w(t)=R(t)$ , we have $${\tilde{R}}_1={\displaystyle {\int }_0^{T/4}{R}_1^2}(t)\text{d}t$$ (16) $${\tilde{R}}_2={\displaystyle {\int }_0^{T/4}{R}_2^2}(t)\mathit{dt}$$ (17)

For simplicity, we can use discrete weighted average $${\tilde{R}}_1=\frac{{\displaystyle {\sum }_{i=1}^N\left\{\left[-A{\omega }_1^2\cos {\omega }_1{t}_i+f(A\cos {\omega }_1{t}_i)\right]w({\omega }_2{t}_i)\right\}}}{{\displaystyle {\sum }_{i=1}^{N}w({\omega }_2{t}_i)}}$$ (18) $${\tilde{R}}_2=\frac{{\displaystyle {\sum }_{i=1}^N\left\{\left[-A{\omega }_2^2\cos {\omega }_2{t}_i+f(A\cos {\omega }_2{t}_i)\right]w({\omega }_2{t}_i)\right\}}}{{\displaystyle {\sum }_{i=1}^{N}w({\omega }_2{t}_i)}}$$ (19)

If we write the trial solution as $$\tilde{u}\text{(}t\text{)}=A\cos \omega t$$ (20)to optimally identify the frequency, this  trequires the residual is minimal $$R(\tilde{u})=-{\omega }^{2}A\cos \omega t+f(\tilde{u})=-{\omega }^2\tilde{u}+f(\tilde{u})\to \min$$ (21)

This implies $$\frac{\mathit{dR}(\tilde{u})}{d\tilde{u}}=-{\omega }^2+\frac{\mathit{df}(\tilde{u})}{d\tilde{u}}=0$$ (22)

We, therefore, obtain another simple frequency formulation, which reads ³⁴ $${\omega }^2=\frac{\mathit{df}}{\mathit{du}}(u=\tilde{u})$$ (23)

In some previous publications, ³⁴ we choose $\tilde{u}=A/2$ $${\omega }^2=\frac{\mathit{df}}{\mathit{dx}}\left(\frac{A}{2}\right)$$ (24)

We can use the weighted average of equation (23) as discussed above.

If we write $$\frac{\mathit{df}}{\mathit{du}}\approx \frac{f(\tilde{u})-f(0)}{\tilde{u}-0}=\frac{f(\tilde{u})}{\tilde{u}}$$ (25)equation (23) becomes ³⁵ $${\omega }^2=\frac{f(\tilde{u})}{\tilde{u}}$$ (26)

This formulation is widely used in the max-min approach ³⁵ $${\omega }^2={\frac{f(A\cos \omega t)}{A\cos \omega t}|}_{\omega t=T/N}$$ (27)

We can also use instead of equation (27) the weighted average of equation (26).

Variational approach

The variational approach^24,29 is to establish a variational principle for equation (1) $$J(u)={\displaystyle {\int }_0^{T/4}\left\{\frac{1}{2}{u^{\prime }}^2-F(u)\right\}\text{d}t}$$ (28)where $\partial F/\partial u=f(u)$ . If we choose an initial guess in the form $$\tilde{u}\text{(}t\text{)}=A\cos \omega t$$ (29)

The functional equation (28) becomes $$J(A)={\displaystyle {\int }_0^{T/4}\left\{\frac{1}{2}{A}^2{\omega }^2{\sin }^2\omega t-F(A\cos \omega t)\right\}\text{d}t}$$ (30)

By setting $$\frac{\mathit{dJ}(A)}{\mathit{dA}}=0$$ (31)we can obtain the needed frequency.

Hamiltonian approach

The Hamiltonian approach ³⁰ is a further extension of the variational approach. In equation (28), ${u^{\prime }}^2/2$ is the kinetic energy and $F(u)$ is the potential energy. According to energy conservation, we have $$\frac{1}{2}{u^{\prime }}^2+F(u)=\frac{1}{2}{u^{\prime }}_0^2+F(u_0)$$ (32)

If we begin with $u\text{(}t\text{)}=A\cos \omega t$ , we have the following equation $$\frac{1}{2}{A}^2{\omega }^2{\sin }^2\omega t+F(A\cos \omega t)-F(A)=0$$ (33)or $${\omega }^2=\frac{2(F(A)-F(A\cos \omega t))}{A^2{\sin }^2\omega t}$$ (34)

Using the weighted average, the frequency can be easily determined. The most simple calculation is to set $\omega t=T/N$ , where T is the period, N is a constant $${\omega }^2={\frac{2(F(A)-F(A\cos \omega t))}{A^2{\sin }^2\omega t}|}_{\omega t=T/N}$$ (35)

Homotopy perturbation method

The homotopy perturbation method^8,9 is to change gradually from a homotopy equation when p = 0 to its original one when p = 1 $$u^{″}+{\omega }^{2}u+p\left[f(u)-{\omega }^{2}u\right]=0$$ (36)

When p = 0, we have a linear oscillator with a frequency of $\omega$ , which is to be determined later. When p = 1, equation (36) becomes the original one.

We can also construct the following homotopy equations in case the Duffing equation’s frequency is known $$u^{″}+u+\epsilon u^3+p\left[f(u)-u-\epsilon u^3\right]=0$$ (37)or $$u^{″}+u+({\epsilon }_0+{\epsilon }_{1}p+{\epsilon }_2{p}^2+\cdots )u^3+p\left[f(u)-u-\epsilon u^3\right]=0$$ (38)where $${\epsilon }_0+{\epsilon }_1+{\epsilon }_2+\cdots =\epsilon$$ (39)

Li and He ¹⁸ suggested a modification of the homotopy perturbation method coupled with the enhanced perturbation method, which is extremely effective for forced nonlinear oscillators; Adamu et al. ²⁰ proposed the parameterized homotopy perturbation method, which gives a high accurate solution for a nonlinear oscillator.

An example

We consider a nonlinear oscillator in the form $$u^{″}+\frac{u^3}{1+u^2}=0\text{,\hspace{1em}}u(0)=A,\hspace{1em}{u}^{\prime }(0)=0$$ (40)

We write equation (40) in the form $$u^{″}+\left(\frac{u^2}{1+u^2}\right)u=0$$ (41)

The coefficient of u in equation (41) is larger than zero, so equation (41) has a periodic solution.

When $u\ll 1$ , equation (40) becomes $$u^{″}+u^3=0,u(0)=A,u^{\prime }(0)=0$$ (42)

The frequency of equation (42) by the homotopy perturbation method is $$\omega =\frac{\sqrt{3}}{2}A$$ (43)

When $u\gg 1$ , equation (41) becomes a linear oscillator $$u^{″}+u=0$$ (44)

Its frequency is $$\omega =1$$ (45)

Now we check the correctness of He’s frequency formulation discussed above. $$f(u)=\frac{u^3}{1+u^2}$$ (46)

It is extremely easy to calculate its derivative with respect to u, that is $$f^{\prime }(u)=\frac{3u^2+u^4}{{(1+u^2)}^2}$$ (47)

Setting u = A/2, we have $${\omega }^2=\frac{\frac{3}{4}{A}^2+\frac{1}{16}{A}^4}{{(1+\frac{1}{4}{A}^2)}^2}$$ (48)

When $A\ll 1$ , equation (48) is equivalent to equation (43); while when $A\gg 1$ , ${\omega }^2=1$ as predicted in equation (45).

The max-min approach leads to the following result $${\omega }^2={\frac{A^2{\cos }^2\omega t}{1+A^2{\cos }^2\omega t}|}_{\omega t=T/12}=\frac{\frac{3}{4}{A}^2}{1+\frac{3}{4}{A}^2}$$ (49)

The variational principle of equation (40) reads $$J(u)={\displaystyle {\int }_0^{T/4}\left\{-\frac{1}{2}{u^{\prime }}^2+\frac{1}{2}{u}^2-\frac{1}{2}\ln (1+u^2)\right\}\text{d}t}$$ (50)

Its Hamiltonian can be written as $$\frac{1}{2}{u^{\prime }}^2+\frac{1}{2}{u}^2-\frac{1}{2}\ln (1+u^2)=\frac{1}{2}{A}^2-\frac{1}{2}\ln (1+A^2)$$ (51)

We choose $u(t)=A\cos \omega t$ and setting $\omega t=T/12$ , we have $$\begin{array}{l}{\omega }^2={\frac{A^2-A^2{\cos }^2\omega t-\ln (1+A^2)+\ln (1+A^2{\cos }^2\omega t)}{A^2{\sin }^2\omega t}|}_{\omega t=T/12}\\ =\frac{A^2-\frac{3}{4}{A}^2-\ln (1+A^2)+\ln (1+\frac{3}{4}{A}^2)}{\frac{1}{4}{A}^2}\\ =\frac{\frac{1}{4}{A}^2-\ln (1+A^2)+\ln (1+\frac{3}{4}{A}^2)}{\frac{1}{4}{A}^2}\end{array}$$ (52)

The accuracy can be improved by weighted average discussed above.

Discussion and conclusions

The above methods are also valid for fractional oscillators.^46–49 We consider a nonlinear oscillator with a fractal derivative^50,51 $$\frac{d^{2}u}{d{t}^{2\alpha }}+·f(u)=0$$ (53)

The fractal derivative $\mathit{du}/d{t}^{\alpha }$ is defined in a fractal time,^50,51 where $\alpha$ is the value of the fractal dimensions. The fractional oscillator implies an intermittent change in time, for example, a stock will not be changed at non-working days, if a week has two non-working days, then the fractal dimensions can be calculated as $$\alpha =\frac{\ln 5}{\ln 7}$$ (54)

For equation (54), we assume that the minimal scale of time is a day, and scales less than one day are not considered. If we observe the stock change in a scale larger than two days, for example, a week or a year, we can consider continuity of a stock change. We assume that the large time scale is s, according to the fractal theory, we have^50,51 $$s=t^{\alpha }$$ (55)

Under a large time scale, equation (53) becomes $$\frac{d^{2}u}{d{s}^2}+·f(u)=0$$ (56)

Equation (56) can be solved analytically by various analytical methods including those discussed above.

The fractal derivative can be used to search for a variational principle for an equation with a damping term $$\frac{d^{2}u}{d{t}^2}+\mu ·\frac{\mathit{du}}{\mathit{dt}}+f(u)=0$$ (57)

The Hamiltonian invariant can be obtained as $$\frac{1}{2}{\left(\frac{\mathit{du}}{\mathit{dt}}\right)}^2+\frac{1}{2}\mu {\left(\frac{\mathit{du}}{d{t}^{1/2}}\right)}^2+F(u)=H_0$$ (58)where $H_0$ is the Hamiltonian invariant, which depends upon the initial conditions and $\partial F/\partial u=f$ .

To conclude, this paper explains some simplest frequency formulae for nonlinear oscillators, the mathematical insight into the formulae reveal all frequency formulae can be improved by choosing a suitable weighted function. This paper can be used as a paradigm for practical applications.