Combination of modified Lindstedt-Poincare and homotopy perturbation methods

Introduction

Mathematical modeling of the physical systems is very important and leads to differential equations many of them are strongly nonlinear such as micro-electromechanical system oscillators, vibration of the trammel pendulum, vibration of a pendulum attached to a rolling wheel, vibration of two-mass spring system etc. Solution of nonlinear differential equations plays an important role in engineering to investigate the physical systems. The methods to solve strongly nonlinear oscillators become inadequate for their analytical treatments. Many researchers have determined the approximate solutions of nonlinear problems by using various methods. Among them the Lindstedt–Poincare (LP) method^1,2 was originally developed to solve

x ¨ + ω 0 2 x + ε f ( x , x , ˙ x ¨ ) = 0 , [ x ( 0 ) = a , x ˙ ( 0 ) = 0 ] , 0 < ε ≪ 1 ,

(1)

Several authors^3–10 extended the method to solve strongly nonlinear problems. Among them Jones ³ presented an approximate technique by introducing a new parameter, α ( ε ) . Such approximate solution is valid for both small and large values of ε . Burton ⁴ and Burton et al. ⁵ presented a modified version of the LP method. Cheung et al. ⁶ further modified the method. However, all approximate solutions obtained by different modified approaches^3–6 are effective for Duffing oscillators with cubical nonlinearity. He^8–10 presented another modification of the LP method. He ¹¹ also presented the homotopy perturbation method for solving equation (1) for large amplitudes of oscillation. In a recent article, ⁷ Alam et al. provided a generalized form of the modified LP ⁶ method and used it to solve various nonlinear oscillators. On the other hand, the homotopy perturbation method has been modified in some recent articles presented by Anjum et al., ¹² He et al., ¹³ Anjum et al., ¹⁴ Ji et al., ¹⁵ and He et al. ¹⁶ Moreover, some approximate techniques are available to solve strongly nonlinear oscillators (Qie et al., ¹⁷ Koochi A and Goharimanesh, ¹⁸ Anjum et al.,^19,20 Rehman et al., ²¹ Suleman et al., ²² Anjum et al., ²³ Li and Nadeem, ²⁴ and Askari et al. ²⁵ ).

The goal of this article is to present a combined technique of the LP and homotopy perturbation methods to solve some strongly nonlinear oscillators. The determination of solutions depends on both parameters α ( ε ) and p → 1 (well-known as a homotopy parameter ¹¹ ). For some nonlinear problems the combined method helps to obtain significantly improved results. An earlier homotopy perturbation method was used as an alternative to the LP method to solve some truly nonlinear oscillators ²⁶ where ω 0 = 0 . A coupling technique of these methods was presented in article ²⁷ . But the present technique is completely different from article ²⁷ . Besides these, the present method is more systematic and easier.

The methods

The modified LP method (Cheung et al. ⁶ )

By a variable transformation τ = ω t , equation (1) becomes

ω 2 x ″ + ω 0 2 x + ε f ( x , ω x ′ , ω 2 x ″ ) = 0 , [ x ( 0 ) = a , x ′ ( 0 ) = 0 ] ,

(2)

where ω is the frequency of this oscillator and primes denote differentiation with respect to τ . Earlier Jones, ³ Burton ⁴ and Burton, and Hamdan ⁵ used the following series for ω

ω 2 = ω 0 2 + ε ω 1 + ε 2 ω 2 + ⋯ ,

(3)

where ω 1 , ω 2 , ⋯ are constants. In article ⁶ , equation (3) was replaced by

ω 2 = ( ω 0 2 + ε ω 1 ) ( 1 + ε 2 ω 2 ω 0 2 + ε ω 1 + ⋯ ) ,

(4)

and then, it was transformed to

ω 2 = ω 0 2 1 − α ( 1 + δ 2 α 2 + ⋯ ) , = ε ω 1 ω 0 2 + ε ω 1 , ε = ω 0 2 α ω 1 ( 1 − α ) ,

(5)

where δ 2 , δ 3 , ⋯ are unknown constants to be determined. Replacing ε , equation (3) was expanded in a series of α as

x = x 0 + α x ∼ 1 + α 2 x ∼ 2 + ⋯ ,

(6)

By substituting equation (5) into equation (2), equation (2) then becomes

( 1 + δ 2 α 2 + ⋯ ) x ″ + ( 1 − α ) x + α f ( x , ω x ′ , ω 2 x ″ ) ω 1 = 0 ,

(7)

Substituting equation (6) into equation (7) and equating coefficients like powers of α , a set of linear differential equations were obtained

x 0 ″ + x 0 = 0 ,

(8)

x ∼ 1 ″ + x ∼ 1 = x 0 − 1 ω 1 f ( x 0 , x 0 ′ , x 0 ″ ) , ⋯ ,

(9)

The initial conditions were replaced by x 0 ( 0 ) = a , x 0 ′ ( 0 ) = 0 , x ∼ 1 ( 0 ) = x ∼ 1 ′ ( 0 ) = x ∼ 2 ( 0 ) = 0 , ⋯ , and x 0 , x ∼ 1 and ω 1 , x ∼ 2 and δ 2 etc. were determined sequentially.

Combination of the modified LP and homotopy perturbation methods

According to the homotopy perturbation technique, equation (7) can be written as

σ x ″ + σ x + p ( − σ x + ( 1 − α ) x + α f ( x , ω x ′ , ω 2 x ″ ) ω 1 ) = 0 ,

(10)

where p is homotopy parameter, ω 2 = ω 0 2 σ ( σ is a dimensionless quantity introduced in article ⁶ ) and

σ = 1 + σ 1 p + σ 2 p 2 + ⋯ , σ 1 = δ 2 α 2 ,

(11)

If p → 1 , equation (10) returns to equation (7). It is noted that the rest of the coefficients σ 2 , σ 3 , ⋯ are not always equal to δ 3 α 3 , δ 4 α 4 , ⋯ for every order of approximations. Therefore, σ becomes slightly different from the result obtained by the modified LP method. ⁶ The solution is also slightly different. A homotopy perturbation solution is considered as

x = x 0 + p x 1 + p 2 x 2 + ⋯ ,

(12)

By substitution of equation (11) and equation (12) into equation (10) and equating coefficients of like powers of p , a set of linear differential equations is obtained. To prevent secular terms, it is assumed that x 1 , x 2 , ⋯ exclude first harmonic terms. Under this assumption all unknown coefficients of two series of x and σ are determined.

Examples (applications)

Some examples are presented in this section to illustrate the accuracy of the proposed method. Numerical results are provided and discussed for some special cases.

Example-1: Duffing oscillator with cubic nonlinearity

The structure for the buckling of a column motion of a particle on a rotating parabola shown in Figure 1 has been considered (Presented in Marinca et al. ²⁸ ). The mass m moves in the horizontal direction only. Neglecting the weight of all but the mass, the equation for the motion of m is

m x ¨ + ( k 1 − 2 p l ) x + ( k 3 − p l 3 ) x 3 + … = 0 ,

(13)

where the spring force is given by

F s p r i n g = k 1 x + k 3 x 3 + … ,

(14)

OPEN IN VIEWER

Figure 1.

Model for the buckling of the column motion of a particle on a rotating parabola.

Now equation (13) can be written in the form

m x ¨ + α 1 x + α 3 x 3 + … = 0 ,

(15)

which is considered as Duffing equation

x ¨ + x + ε x 3 = 0 ,

(16)

Let us consider the initial condition as x ( 0 ) = a , x ˙ ( 0 ) = 0 .

For equation (16), f = x 3 and equation (10) becomes

x ″ + x + p ( − x + ( ( 1 − α ) x + α x 3 ω 1 ) / σ ) = 0 ,

(17)

According to the processes discussed in the previous section, the linear differential equations for x 1 , x 2 , ⋯ have been obtained

x 0 ″ + x 0 = 0 ,

(18)

x 1 ″ + x 1 = α ( − x 0 + x 0 3 ω 1 ) ⋯ ,

(19)

x 2 ″ + x 2 = α ( − x 1 + 3 x 0 2 x 1 ω 1 ) − σ 1 ( ( 1 − α ) x 0 + x 0 3 ω 1 ) , ⋯ ,

(20)

Solving the above set of equations with the initial conditions described in the previous section, we obtained the following results

σ 1 = − α 2 24 , σ 2 = − α 4 576 , σ 3 = 7 α 4 13824 − α 6 6912 , σ 4 = 7 α 6 82944 − 5 α 8 331776 , ⋯ ,

(21)

and

x 0 = a cos τ , x 1 = a α 24 ( cos 3 τ − cos τ ) ,

x 2 = a α 3 72 ( cos 3 τ − cos τ ) + a α 2 24 ( cos 5 τ − cos τ ) ,

(22)

x 3 = a ( 3 − α 2 ) 6912 ( cos 3 τ − cos τ ) − a α 4 6912 ( cos 5 τ − cos τ ) − a α 3 13824 ( cos 7 τ − cos τ ) , ⋯ ,

Substituting p = 1 together with the values of σ 1 , σ 2 , ⋯ , equation (22) becomes

σ = 1 + ∑ i = 1 4 σ i = 1 − α 2 24 − 17 α 4 13824 − 5 α 6 82944 − 5 α 8 331776 .

(23)

By the modified LP method, ⁶ the following results was obtained

δ 2 = − α 2 24 , δ 3 = 0 , δ 4 = − 17 α 4 13824 , δ 5 = 0 , ⋯ ,

σ = 1 + ∑ i = 1 4 δ i + 1 = 1 − α 2 24 − 17 α 4 13824 .

(24)

In both relations of σ , the results are the same until α 4 and it is the correct result. If the next terms σ 5 and δ 6 are added respectively in the series equation (23) and Eq. (24), the results will be the same until α 6 and later series (presented by equation (24)) do not contain more than α 6 ordered terms. On the contrary, the former contains more two terms up to α 10 . Again, all the coefficients are negative. It is observed that all approximate results of σ are smaller than its exact value and the former (obtained by a combined technique) provides better results for this oscillator.

Example-2: nonlinear oscillator with inertia-type cubic nonlinearity

The motion of a ring of mass m sliding freely on a wire described by the parabola y = q x 2 which rotates with a constant angular velocity Ω about the y-axis as shown in Figure 2 (Nayfeh and Mook ² ). The equation describing the motion of the ring is (Nayfeh and Mook, ² Marinca and Herinasu ²⁸ ).

x ¨ + ω 2 x = − 4 q x ( x x ¨ + x ˙ 2 ) ,

(25)

where ω 2 = 2 g q − Ω 2 and the initial conditions are x ( 0 ) = A , x ˙ ( 0 ) = 0 . OPEN IN VIEWER

Figure 2.

Motion of a particle on a rotating parabola.

Equation (25) can be written as ²⁹

x ¨ + x + ε ( x 2 x ¨ + x ˙ 2 x ) = 0 ,

(26)

with its transformed form

ω 2 x ″ + x + ε ( ω 2 x 2 x ″ + ω 2 x ′ 2 x ) = 0 ,

(27)

or,

x ″ + ω − 2 x + ε ( x 2 x ″ + x ′ 2 x ) = 0 ,

(28)

For this oscillator, by choosing the following series (see article ⁷ for details)

ω − 2 = 1 + ε ω 1 + ε 2 ω 2 + ⋯ ,

(29)

and the corresponding equation of equation (7) becomes

( 1 − α ) x ″ + σ x + α ω 1 ( x 2 x ″ + x ′ 2 x ) = 0 ,

(30)

where σ follows equation (11).

Now adding a homotopy parameter it becomes

σ x ″ + σ x + p ( − σ x ″ + ( 1 − α ) x ″ + α ω 1 ( x 2 x ″ + x ′ 2 x ) ) = 0 ,

(31)

For equation (31), the following results have been obtained

σ 1 = − α 2 8 , σ 2 = − α 4 64 , σ 3 = − 5 α 4 512 − α 6 256 , σ 4 = − 5 α 4 512 − 5 α 8 4096 , ⋯ ,

and then,

σ = 1 + ∑ i = 1 4 σ i = 1 − α 2 8 − 13 α 4 512 − 9 α 6 1024 − 5 α 8 4096 ,

(32)

On the other hand, the calculated result by the modified LP method ⁷ is

σ = 1 + ∑ i = 1 4 δ i + 1 = 1 − α 2 8 − 13 α 4 512 ,

(33)

All these results are similar to those obtained for the cubic Duffing oscillator (i.e., the coefficients of α 3 , α 5 , ⋯ vanish). For both, oscillators δ 2 i + 1 = 0 , i = 0 , 1 , 2 , ⋯ and coefficients are much smaller than 1 .

Example-3: duffing oscillator with quintic nonlinearity

Let us consider a highly nonlinear oscillator, that is, Duffing oscillator with quintic nonlinearity

x ¨ + ω 0 2 x + ε x 5 = 0 ,

(34)

For this oscillator the following results have been obtained

σ 1 = − 13 α 2 120 , σ 2 = 2 α 3 225 − 169 α 4 14400 , σ 3 = 1673 α 4 575000 + 13 α 3 4500 − 2197 α 6 864000 ,

σ 4 = − 253 α 5 270000 + 11411 α 6 10368000 + 169 α 7 162000 − 28561 α 8 41472000 , ⋯ ,

and

σ = 1 + ∑ i = 1 4 σ i = 1 − 13 α 2 120 + 2 α 3 225 − 5087 α 4 576000 + 527 α 5 270000 − 14953 α 6 10368000 + 169 α 7 162000 − 28561 α 8 41472000 ,

(35)

On the other hand, the calculated result by the modified LP method ⁷ is

σ = 1 + ∑ i = 1 4 δ i + 1 = 1 − 13 α 2 120 + 2 α 3 225 − 5087 α 4 576000 + 527 α 5 270000 ,

(36)

The coefficients of equation (36) are also much smaller than 1 , but larger than those of equation (23) and equation (24). Moreover, all δ 2 i + 1 , i = 0 , 1 , 2 , ⋯ are positive. For this reason, the results oscillate and gradually approach its exact value. However, the errors of results obtained by equation (35) are much smaller than those obtained by equation (36).

Example-4: duffing oscillator with septic nonlinearity

Let us consider another highly nonlinear oscillator, that is, Duffing oscillator with septic nonlinearity

x ¨ + ω 0 2 x + ε x 7 = 0 ,

(37)

For this oscillator the following results have been obtained

σ 1 = − 99 α 2 560 , σ 2 = 141 α 3 4480 − 9801 α 4 313600 , σ 3 = 695357 α 4 175616000 + 41877 α 5 2508800 − 970299 α 6 87808000 ,

σ 4 = − 2917107 α 5 614656000 + 20131893 α 6 24586240000 + 1381941 α 7 140492800 − 96059601 α 8 19668992000 , ⋯ ,

and

σ = 1 + ∑ i = 1 4 σ i = 1 − 99 α 2 560 + 141 α 3 4480 − 4793203 α 4 175616000 + 3671379 α 5 307328000 − 251551827 α 6 24586240000 + 1381941 α 7 140492800 − 96059601 α 8 19668992000 ,

(38)

On the other hand, the calculated result by the modified LP method ⁷ is

σ = 1 + ∑ i = 1 4 δ i + 1 = 1 − 99 α 2 560 + 141 α 3 4480 − 4793203 α 4 175616000 ,

(39)

The coefficients of equation (39) are also much smaller than 1 but larger than those of equation (23) and equation (24). Moreover, all δ 2 i + 1 , i = 0 , 1 , 2 , ⋯ are positive. For this reason, the results oscillate and gradually approach its exact value. However, the errors of results obtained by equation (38) are much smaller than those obtained by equation (39).

Convergency of the α -series

For a lot of nonlinear oscillators, the α -series of the modified LP method ⁷ converge rapidly when α ≪ 1 (i.e., in the case of small oscillations). In the case of large amplitude of oscillations (i.e., α → 1 ) the series also converges since the coefficients of the series, δ i , i = 2 , 3 , ⋯ gradually decrease. For some of the oscillators, the coefficients δ i , i = 3 , 5 , ⋯ vanish. According to the recently presented homotopy perturbation method, ¹⁶ the coefficients ω i , i = 3 , 5 , ⋯ vanish for oscillator presented by equation (16). In these situations, the frequency of oscillation is same for second and third approximations; and then for fourth and fifth approximations. On the other hand, the coefficients δ i , i = 3 , 5 , ⋯ are different in sign of the coefficients δ i , i = 2 , 4 , ⋯ for oscillators presented by equation (34) and equation (37). In these situations, approximate results of frequencies oscillate. For these nonlinear oscillators, the calculated frequencies obtained by the homotopy perturbation method ¹⁶ also oscillates. Moreover, homotopy perturbation method ¹⁶ fails to measure frequency as well as solution of oscillator presented by equation (26) especially in the case of amplitude of oscillation. All these problems have been satisfactorily removed in this article.

Results and discussion

Recently the modified Lindstedt–Poincare method ⁶ has been generalized for solving strongly nonlinear oscillators. ⁷ The generalized method covers wide varieties of nonlinear problems where ω 0 ≠ 0 . (It is noted that the modified technique ⁶ is useless for a class of nonlinear oscillators, that is, second example equation (26)). Moreover, the generalized version ⁷ has been combined with the homotopy perturbation technique for some truly nonlinear oscillators where ω 0 = 0 . In this article the method has been combined with homotopy perturbation for all values of ω 0 whether it vanishes or not. Earlier truly nonlinear oscillators were solved by the homotopy method, ¹¹ but with a different approach. The new technique is more systematic than all the past methods mentioned above. For large amplitudes of oscillations α → 1 and errors of frequency become the largest. To compare the results with existing results, the ratios of frequency with their exact frequency are presented below.

For the cubic Duffing oscillator equation (16), the ratios respectively become:

By the present method with equation (23): lim ε a 2 → ∞ . ω ω E x a c t ≅ 1.0000008 .

and by the modified LP method ⁶ with equation (24): . lim ε a 2 → ∞ ω ω E x a c t ≅ 1.00004 .

For the nonlinear oscillators equation (26), the ratios respectively become

by the present method with equation (32): . lim ε a 2 → ∞ ω ω E x a c t ≅ 0.98256 ,

and by the modified LP method, ⁷ with equation (33): . lim ε a 2 → ∞ ω ω E x a c t ≅ 0.97676 ,

For the highly nonlinear oscillators equation (34), the ratios respectively become

by the present method with equation (35): lim ε a 4 → ∞ . ω ω E x a c t ≅ 1.000095 ⋯ ,

and by the modified LP method, ⁷ with equation (36) : lim ε a 4 → ∞ ω ω E x a c t ≅ 1.000704 ⋯ .

For the Duffing oscillator with cubic nonlinearity, the second, third, fourth, and fifth approximated frequencies with their percentage error are given in Table 1.

OPEN IN VIEWER

Table 1.

Comparison of the frequencies among the present method, modified Lindstedt–Poincare method, ⁷ and the homotopy method ¹⁶ for the Duffing oscillator with cubic nonlinearity for a = 100 .

Approx.	Present method (% Err.)	MLP method 7 (% Err.)	Homotopy method 16 (%Err.)
2nd	84.7853 (0.06820)	84.7853 (0.06820)	84.7013 (0.03087)
3rd	84.7085 (0.02243)	84.7853 (0.06820)	84.7013 (0.03087)
4th	84.7245 (0.00354)	84.7309 (0.00401)	84.7281 (0.00073)
5th	84.7275 (0.00008)	84.7309 (0.00401)	84.7281 (0.00073)

Moreover, for the highly nonlinear oscillators equation (34), the second, third, fourth, and fifth approximated frequencies with their percentage error are given in Table 2.

OPEN IN VIEWER

Table 2.

Comparison of the frequencies among the present method, modified Lindstedt–Poincare method, ⁷ and the homotopy method ¹⁶ for the quintic oscillator ( x ¨ + ω 0 2 x + ε x 5 = 0 ) for a = 100 .

Approx.	Present method (% Err.)	MLP method 7 (% Err.)	Homotopy method 16 (% Err.)
2nd	7465.2 (0.0420)	7465.2 (0.0420)	7400.95 (0.90234)
3rd	7465.44 (0.03883)	7502.31 (0.4548)	7455.8541 (0.16718)
4th	7467.57 (0.01031)	7465.44 (0.0388)	7474.9977 (0.089146
5th	7469.05 (0.0095)	7473.60 (0.0704)	7468.15 (0.00254)

The present solution along with numerical, Modified Lindstedt–Poincare ⁶ and generalized Lindstedt–Poincare ⁷ solutions are presented in Figures 3–7. Also, the frequency analysis of the present method with existing methods is provided in Figures 8 and 9.OPEN IN VIEWER

OPEN IN VIEWER

Figure 4.

Comparison between the time versus displacement results obtained by the present method, numerical method, and Cheung et al. ⁶ method for the Duffing oscillator with cubic nonlinearity when a = 100 , ε = 1 .

OPEN IN VIEWER

Figure 5.

Comparison between the time versus displacement results obtained by the present method, numerical method, Cheung et al. ⁶ method, and Alam et al. ⁷ method for equation (26) when a = 0.7 , ε = 1 .

OPEN IN VIEWER

Figure 6.

Comparison between the time versus displacement results obtained by the present method, numerical method, and Alam et al. ⁷ method for equation (26) when a = 5.0 , ε = 1 .

OPEN IN VIEWER

Figure 7.

Comparison between the time versus displacement results obtained by the present method, numerical method, Cheung et al. ⁶ method, and Alam et al. ⁷ method for the quintic oscillator ( x ¨ + ω 0 2 x + ε x 5 = 0 ) when a = 0.7 , ε = 1 .

OPEN IN VIEWER

Figure 8.

Comparison of calculated frequencies measured by the present technique and earlier modified technique (Alam et al. ⁷ ), homotopy method ¹⁶ together with the numerical results for the quintic oscillator ( x ¨ + ω 0 2 x + ε x 5 = 0 ) when a = 100 , ε = 1 .

OPEN IN VIEWER

Figure 9.

Comparison of calculated frequencies measured by the present technique and earlier modified technique (Alam et al. ⁷ ) together with the numerical results for the septic oscillator ( x ¨ + ω 0 2 x + ε x 7 = 0 ) when a = 10.0 , ε = 1 .

All these figures show that the present results have good agreement with the numerical results; on the other hand, in some cases, the present results are better (Figure 3, Figures 5–9) than other existing results. Also, for a large amplitude, the present results are converging but other results are diverging, that is, the results obtained by Cheung et al. ⁶ diverge for equation (26) and the results obtained by Alam et al. ⁷ diverge for Duffing cubic, quintic, and septic oscillators.

Conclusion

Earlier modified LP method and homotopy perturbation method were presented for solving strongly nonlinear oscillators especially in case of large amplitude of oscillation. Based on these methods, a combined technique has been presented which provides more satisfactory results. The necessity of present technique has been illustrated by considering four different nonlinear oscillators. Calculated frequencies by the modified LP method of cubical Duffing oscillator (i.e., equation (24)) gradually closes to numerical result with increasing order of approximation. But present frequency (equation (23)) closes faster than that of the modified LP method and a recent homotopy perturbation method. ¹⁶ For the quintic oscillator, coefficients δ i , i = 3 , 5 , ⋯ are different in sign of δ i , i = 2 , 4 , ⋯ . It is noted that δ i = 0 , i = 3 , 5 , ⋯ for cubical oscillator. So, major problem arises in case of quintic oscillator. In this situation error of the frequency (equation (36), which is obtained by modified LP method) oscillates. Similar problem arises in septic oscillator. Homotopy perturbation method ¹⁶ measures also oscillating results for these oscillators. But present method provides more accurate result. Moreover, the present method covers various nonlinear oscillators involving x , x ˙ and x ¨ (e.g., equation (26)). It is noted that the homotopy perturbation method ¹⁶ is useless for solving equation (26). Homotopy perturbation method ¹⁶ has another limitation. According to this method, n-th ( n = 2 , 3 , ⋯ ) approximate frequency are calculated from n-th order algebraic equation. On the contrary, modified LP method and present method are straightforward. However, the present method has a limitation for solving the quadratic oscillators. It measures similar results of modified LP method. ⁶