Analytical periodic solution for solving nonlinear vibration equations

1. Introduction

Finding an analytical solution for nonlinear vibration problems has been a topic of intensive research for many years.

In the past recent decades some powerful and new analytical methods have been considered by many researchers, such as: homotopy perturbation method [1], energy balance method [2-5], amplitude frequency formulation [6], iteration perturbation method [7], max-min approach [8-10], variational approach method [11-16], homotopy analysis method [17], and the other analytical and numerical methods [18-24].

Among these methods, Variational Approach Method (VAM) is considered to solve the nonlinear vibration equations in this paper.

The paper has been collocated as follows.

The basic idea of He’s variational approach is introduced in section 2. In section 3 for numerical solution of nonlinear equation the algorithm of Runge-Kutta method is described. Then in section 4 applications of He’s variational approach method have been studied, to demonstrate the applicability and preciseness of the method. In section 5 some comparisons between analytical and numerical solutions are presented. Eventually we show that VAM can converge to a precise cyclic solution for nonlinear systems.

2. Basic concept of variational approach method

He suggested a variational approach which is different from the known variational methods in open literature [14]. Hereby we give a brief introduction of the method:

Its variational principle can be easily established utilizing the semi-inverse method:

where $T$ is period of the nonlinear oscillator, $\raisebox{1ex}{$\partial F$}\!\left/ \!\raisebox{-1ex}{$\partial u$}\right.=f$. Assume that its solution can be expressed as:

where $A$ and $\omega$ are the amplitude and frequency of the oscillator, respectively. Substituting Eq. (3) into Eq. (2) results in:

Applying the Ritz method, He requires:

But with a careful inspection, for most cases He found that:

Thus He modified conditions Eq. (5) and Eq. (6) into a simpler form:

from which the relationship between the amplitude and frequency of the oscillator can be obtained.

3. Basic concept of Runge-Kutta method

For the numerical approach to verify the analytic solution, the fourth order RK (Runge-Kutta) method has been used. This iterative algorithm is written in the form of the following formulae for the second-order differential equation:

where $\Delta t$ is the increment of the time and $h_1$, $h_2$, $h_3$ and $h_4$ are determined from the following formulae:

The numerical solution starts from the boundary at the initial time, where the first value of the displacement function and its first-order derivative are determined from initial condition. Then with a small time increment $\Delta t$ the displacement function and its first-order derivative at the new position can be obtained using Eq. (10). This process continues to the end of the time limit.

4. Application

In order to show the accuracy of the new approach method, we will apply the method for the following two examples.

4.1. Example 1

We consider the physical model of nonlinear equation in the following figure with $F\left(t\right)=F_0\mathrm{s}\mathrm{i}\mathrm{n}{\omega }_{0}t$ indicated in Figure 1.

The motion equation is:

11

$\ddot{\theta }+\frac{4k}{3m}\mathrm{s}\mathrm{i}\mathrm{n}\theta -\frac{3F_0}{ml}\mathrm{s}\mathrm{i}\mathrm{n}{\omega }_{0}t=0,\theta \left(0\right)=A,\dot{\theta }\left(0\right)=0.$

Fig. 1The physical model of nonlinear equation

This equation is known as Mathieu equation or the system with dependent coefficients to time, in which $\theta$ and $t$ are generalized dimensionless displacements and time variables, respectively.

The approximation $\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)=\theta -\left(1/6\right){\theta }^3+\left(1/120\right){\theta }^5$ is used.

Its variational formulation can be readily obtained from Eq. (11) as follows:

12

$J\left(\theta \right)={\int }_0^t\left(\frac{1}{2}{\ddot{\theta }}^2+\frac{2}{3}\frac{k}{m}{\theta }^2-\frac{1}{18}\frac{k}{m}{\theta }^4+\frac{1}{540}\frac{k}{m}{\theta }^6-\frac{3F_0\sin \left({\omega }_{0}t\right)}{ml}\theta \right)dt.$

Choosing the trial function $\theta \left(t\right)=A\mathrm{c}\mathrm{o}\mathrm{s}\left(\omega t\right)$ into Eq. (12) we obtain:

13

$J\left(A\right)={\int }_0^{T/4}\left(\begin{array}{c}\frac{1}{2}{A}^2{\omega }^2{\mathrm{s}\mathrm{i}\mathrm{n}}^2\left(\omega t\right)+\frac{2}{3}\frac{k}{m}{A}^2{\mathrm{c}\mathrm{o}\mathrm{s}}^2\left(\omega t\right)-\frac{1}{18}\frac{k}{m}{A}^4{\mathrm{c}\mathrm{o}\mathrm{s}}^4\left(\omega t\right)\\ +\frac{1}{540}\frac{k}{m}{A}^6{\mathrm{c}\mathrm{o}\mathrm{s}}^6\left(\omega t\right)-\frac{3F_0\sin \left({\omega }_{0}t\right)}{ml}Acos\left(\omega t\right)\end{array}\right)dt.$

The stationary condition with respect to $A$ leads to:

14

$\frac{\partial J}{\partial A}={\int }_0^{T/4}\left(\begin{array}{c}A{\omega }^2{\mathrm{s}\mathrm{i}\mathrm{n}}^2\left(\omega t\right)+\frac{4}{3}\frac{k}{m}A{\mathrm{c}\mathrm{o}\mathrm{s}}^2\left(\omega t\right)-\frac{2}{9}\frac{k}{m}{A}^3{\mathrm{c}\mathrm{o}\mathrm{s}}^4\left(\omega t\right)\\ +\frac{1}{90}\frac{k}{m}{A}^5{\mathrm{c}\mathrm{o}\mathrm{s}}^6\left(\omega t\right)-\frac{3F_0\sin \left({\omega }_{0}t\right)}{ml}\cos \left(\omega t\right)\end{array}\right)dt=0,$

or:

15

$\frac{\partial J}{\partial A}={\int }_0^{\pi /2}\left(\begin{array}{c}A{\omega }^2{\mathrm{s}\mathrm{i}\mathrm{n}}^{2}t+\frac{4}{3}\frac{k}{m}A{\mathrm{c}\mathrm{o}\mathrm{s}}^{2}t-\frac{2}{9}\frac{k}{m}{A}^3{\mathrm{c}\mathrm{o}\mathrm{s}}^{4}t\\ +\frac{1}{90}\frac{k}{m}{A}^5{\mathrm{c}\mathrm{o}\mathrm{s}}^{6}t-\frac{3F_0\sin \left({\omega }_{0}t\right)}{ml}cost\end{array}\right)dt=0.$

Solving Eq. (15) according to $\omega$ we have:

16

$\omega {}^2=\frac{{\int }_0^{\frac{\pi }{2}}\left(\frac{4}{3}\frac{k}{m}A{\mathrm{c}\mathrm{o}\mathrm{s}}^{2}t-\frac{2}{9}\frac{k}{m}{A}^3{\mathrm{c}\mathrm{o}\mathrm{s}}^{4}t+\frac{1}{90}\frac{k}{m}{A}^5{\mathrm{c}\mathrm{o}\mathrm{s}}^{6}t-\frac{3F_0\sin \left({\omega }_{0}t\right)}{ml}\mathrm{c}\mathrm{o}\mathrm{s}t\right)dt}{{\int }_0^{\frac{\pi }{2}}A{\mathrm{s}\mathrm{i}\mathrm{n}}^{2}tdt}.$

Then we have:

17

$\omega {}_{VAM}=\frac{1}{12}\sqrt{\frac{1728F_0\sin \left(\frac{1}{2}\pi {\omega }_0\right)-1728F_0{\omega }_0+kAl\pi \left({\omega }_0^2-1\right)\left(192+A^4-24A^2\right)}{\left(m{\omega }_0^2-m\right)lA\pi }}.$

According to Eqs. (3) and (17) we can obtain the following approximate solution:

18

$\theta \left(t\right)=A\cos \left(\frac{1}{12}\sqrt{\frac{1728F_0\sin \left(\frac{1}{2}\pi {\omega }_0\right)-1728F_0{\omega }_0+kAl\pi \left({{\omega }_0}^2-1\right)\left(192+A^4-24A^2\right)}{\left(m{{\omega }_0}^2-m\right)lA\pi }}t\right).$

In Figures 3 and 4 we compared variational approach method with numerical solution for different parameters.

4.2. Example 2

We consider the motion equation of the pendulum with harmonic stringer point in Figure 2.

The solution of Mathieu equation of pendulum with harmonic stringer point according to the variational approach method is:

19

$\ddot{\theta }+\left(\frac{g}{l}-\frac{{\omega }_0^{2}Y}{l}\mathrm{c}\mathrm{o}\mathrm{s}{\omega }_{0}t\right)\mathrm{s}\mathrm{i}\mathrm{n}\theta =0,\mathrm{}\mathrm{}\theta \left(0\right)=A,\dot{\theta }\left(0\right)=0.$

This equation is known as Mathieu equation or the system with dependent coefficients to time, in which $\theta$ and $t$ are generalized dimensionless displacements and time variables, respectively.

The approximation $\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)=\theta -\left(1/6\right){\theta }^3$ is used.

Its variational formulation can be readily obtained from Eq. (19) as follows:

20

$J\left(\theta \right)={\int }_0^t\left(\frac{1}{2}{\ddot{\theta }}^2+\frac{1}{2}\left(\frac{g}{l}-\frac{{\omega }_0^{2}Y\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{0}t\right)}{l}\right){\theta }^2+\frac{1}{24}\left(\frac{g}{l}-\frac{{\omega }_0^{2}Y\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{0}t\right)}{l}\right){\theta }^4\right)dt.$

Fig. 2Pendulum with harmonic stringer point: y(t)=Ycosω0 t

Choosing the trial function $\theta \left(t\right)=A\mathrm{c}\mathrm{o}\mathrm{s}\left(\omega t\right)$ into Eq. (18) we obtain:

21

$J\left(A\right)={\int }_0^{T/4}\left(\begin{array}{c}\frac{1}{2}{A}^2{\omega }^2{\mathrm{s}\mathrm{i}\mathrm{n}}^2\left(\omega t\right)+\frac{1}{2}{A}^2\left(\frac{g}{l}-\frac{{{\omega }_0}^{2}Y\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{0}t\right)}{l}\right){\mathrm{c}\mathrm{o}\mathrm{s}}^2\left(\omega t\right)\\ +\frac{1}{24}{A}^4\left(\frac{g}{l}-\frac{{{\omega }_0}^{2}Y\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{0}t\right)}{l}\right){\mathrm{c}\mathrm{o}\mathrm{s}}^4\left(\omega t\right)\end{array}\right)dt.$

The stationary condition with respect to $A$ leads to:

22

$\frac{\partial J}{\partial A}={\int }_0^{T/4}\left(\begin{array}{c}A{\omega }^2{\mathrm{s}\mathrm{i}\mathrm{n}}^2\left(\omega t\right)+A\left(\frac{g}{l}-\frac{{\omega }_0^{2}Y\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{0}t\right)}{l}\right){\mathrm{c}\mathrm{o}\mathrm{s}}^2\left(\omega t\right)\\ +\frac{1}{6}A\left(\frac{g}{l}-\frac{{\omega }_0^{2}Y\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{0}t\right)}{l}\right){\mathrm{c}\mathrm{o}\mathrm{s}}^4\left(\omega t\right)\end{array}\right)dt=0,$

or:

23

$\frac{\partial J}{\partial A}={\int }_0^{\pi /2}\left(\begin{array}{c}A{\omega }^2{\mathrm{s}\mathrm{i}\mathrm{n}}^{2}t+A\left(\frac{g}{l}-\frac{{\omega }_0^{2}Y\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{0}t\right)}{l}\right){\mathrm{c}\mathrm{o}\mathrm{s}}^{2}t\\ +\frac{1}{6}A\left(\frac{g}{l}-\frac{{\omega }_0^{2}Y\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{0}t\right)}{l}\right){\mathrm{c}\mathrm{o}\mathrm{s}}^{4}t\end{array}\right)dt=0.$

Solving Eq. (23) according to $\omega$ we have:

24

$\omega {}^2=\frac{{\int }_0^{\frac{\pi }{2}}\left(\left(\frac{g}{l}-\frac{{\omega }_0^{2}Y\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{0}t\right)}{l}\right){\mathrm{c}\mathrm{o}\mathrm{s}}^{2}t+\frac{1}{6}A\left(\frac{g}{l}-\frac{{\omega }_0^{2}Y\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{0}t\right)}{l}\right){\mathrm{c}\mathrm{o}\mathrm{s}}^{4}t\right)dt}{{\int }_0^{\frac{\pi }{2}}A{\mathrm{s}\mathrm{i}\mathrm{n}}^{2}tdt}.$

Fig. 3Comparison of analytical solution of θ based on time with the numerical solution for cases: a) L=0.5 m, m=20 kg, k=800 N/m, F0=1 N, ω0=2 rad/sec, A=π/6, b) L=1 m, m=10 kg, k=1200 N/m, F0=1 N, ω0=2 rad/sec, A=π/3

a)

b)

Fig. 4Comparison of phase plan of the analytical solution with the numerical solution for cases: a) L=0.5 m, m=20 kg, k=800 N/m, F0=1 N, ω0=2 rad/sec, A=π/6, b) L=1 m, m=10 kg, k=1200 N/m, F0=1 N, ω0=2 rad/sec, A=π/3

a)

b)

Then we have:

25

$\omega {}_{VAM}=\frac{\sqrt{2}}{4}\sqrt{\frac{128Y\left(A^2+\frac{1}{2}{\omega }_0^2-8\right){\omega }_0\sin \left(\frac{1}{2}\pi {\omega }_0\right)-g\pi \left({\omega }_0^2-16\right)\left({\omega }_0^2-4\right)\left(A^2-8\right)}{l\left({\omega }_0^4-20{\omega }_0^2+64\right)\pi }}.$

According to Eqs. (3) and (25) we can obtain the following approximate solution:

26

$\theta \left(t\right)=A\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\sqrt{2}}{4}\sqrt{\frac{128Y\left(A^2+\frac{1}{2}{\omega }_0^2-8\right){\omega }_0\sin \left(\frac{1}{2}\pi {\omega }_0\right)-g\pi \left({\omega }_0^2-16\right)\left({\omega }_0^2-4\right)\left(A^2-8\right)}{l\left({\omega }_0^4-20{\omega }_0^2+64\right)\pi }}t\right).$

For this example we consider different parameters to show the accuracy of the method.

Figures 5 and 6 show the excellent agreement of the variational approach method with the numerical one.

Fig. 5Comparison of analytical solution of θ t based on time with the numerical solution: a) L=0.5 m, ω0=1rad/sec, Y=0.25 m, g=9.81m/s2, A=π/6, b) L=1 m, ω0=1 rad/sec, Y=0.5 m, g=9.81 m/s2, A=π/4

a)

b)

Fig. 6Comparison of phase plane of the analytical solution with the numerical solution for: a) L=0.5 m, ω0=1 rad/sec, Y=0.25 m, g=9.81 m/s2, A=π/6, b) L=1 m, ω0=1 rad/sec, Y=0.5 m, g=9.81 m/s2, A=π/4

a)

b)

5. Results and discussions

Some diagrams are illustrated to show the accuracy of the variational approach method with the numerical solution using Runge-Kutta method. Figure 3 represents a comparison of analytical solution of $\theta \left(t\right)$ based on time with the numerical solution for two cases (a): $L=0.5\mathrm{}\mathrm{m}$, $m=20\mathrm{k}\mathrm{g}$, $k=800\mathrm{}\mathrm{N}/\mathrm{m}$, $F_0=1\mathrm{N}$, ${\omega }_0=2\mathrm{}\mathrm{r}\mathrm{a}\mathrm{d}/\mathrm{s}\mathrm{e}\mathrm{c}$, $A=\pi /\text{6}$, (b): $L=1\mathrm{}\mathrm{m}$, $m=10\mathrm{}\mathrm{k}\mathrm{g}$, $k=1200\mathrm{}\mathrm{N}/\mathrm{m}$, $F_0=1\mathrm{N}$, ${\omega }_0=2\mathrm{}\mathrm{r}\mathrm{a}\mathrm{d}/\mathrm{s}\mathrm{e}\mathrm{c}$, $A=\pi /\text{3}$. Figure 4 shows the comparisons of analytical solution of $\dot{\theta }$ based on $\theta \left(t\right)$ with the numerical solution again for two important cases (i) and (ii). In example 2, Figure 5 represents the comparison of analytical solution of $\theta \left(t\right)$ based on time with the numerical solution for two cases (a) and (b). It shows that the behavior of the oscillation is periodic. And Figure 6 represents the comparison of $\dot{\theta }$ versus $\theta \left(t\right)$ of the analytical solution with the numerical solution for different parameters in two cases to show the accuracy of the method. It has been indicated that VAM has excellent agreement with the numerical solution. It is a simple method and easy to apply to any kind of nonlinear vibration problems.

6. Conclusions

In this study a new mathematical tool has been successfully used to study nonlinear vibration equations. We considered two highly nonlinear problems in this paper, the results obtained from the VAM and Runge-Kutta method have an excellent agreement. It has been shown that the obtained solutions are valid for the whole domain. In VAM just one iteration leads us to a high accuracy of the solution. It has been shown that the VAM can be easily extended to any kind of nonlinear vibration equations.