Sound & Vibration

DOI: 10.32604/sv.2022.014166

ARTICLE

Damped Mathieu Equation with a Modulation Property of the Homotopy Perturbation Method

Yusry O. El-Dib and Nasser S. Elgazery*

Department of Mathematics, Faculty of Education, Ain Shams University, Roxy, Cairo, Egypt
*Corresponding Author: Nasser S. Elgazery. Email: nasser522000@gmail.com
Received: 06 September 2020; Accepted: 25 June 2021

1  Introduction

In wide engineering and physical applications, the nonlinear oscillators exist. Also, the parametric excitation takes place when a modifying physical parameter, such as a moment of stiffness or inertia, acts in a forcing model. This excitement yields a variable time coefficient, commonly an oscillation, in the governing system of motion. On the other hand, an external excitation outcome acting as an inhomogeneous part in the model of motion. Furthermore, minor parametric excitement produces a major response when the frequency of the excitement is far from the fundamental resonance, as shown in [1–5].

A classical example of parametric excitation is the swinging pendulum with oscillating support. The equation of motion describing the model is the well-known Mathieu equation. In 1868 Mathieu studied the vibration of elliptical membranes [6]. Consequently, he introduced the Mathieu equation that is an example of a linear differential equation (LDE) with parametric excitation. The Mathieu equation has application to the dynamics of passive towed arrays in submarines, as well as serving as a useful model for many interesting problems in physics, biology, applied mathematics, and engineering mechanics fields [7]. For some of the non-linear variations of the Mathieu, the equation has been presented in [8,9]. Moreover, the oscillations of the mechanical systems under the action of an oscillatory external force may reveal a Duffing problem, for instance, see references [10–14]. Recently, Moatimid [15] attempted to study the stability analysis of a parametric Duffing oscillator. In this investigation Moatimid showed that the cubic stiffness parameter and the damped parameter have a destabilizing influence, however, the parametric and natural frequencies are of stabilizing influences.

The main target in the present work is how to achieve accurate approximate solutions of the nonlinear oscillators with highly strong nonlinearity. In recent centuries, many analytical approaches were developed to work out the periodic motion of nonlinear oscillators, such as the averaging method, perturbation methods, harmonic balance method, and the generalized harmonic method. The classical perturbation procedure depends on small parameters and chooses unsuitable small parameters that can lead to wrong solutions [16]. Therefore, a new perturbation technique was first proposed by He et al. [17–20]. This technique is named as the HPM, which represents a combination of the Homotopy analysis and classical perturbation methods. It has a full promise of the traditional perturbation techniques. The major property of the HPM is in its ability and flexibility to deal with many types of linear and nonlinear differential equations conveniently and accurately. Further, the HPM provides us with an appropriate direction to calculate an approximate or an analytic solution to several models arising in different fields. He [21] was built the most two considerable steps in the criteria of the HPM with a suitable initial guess and suggested an alternative approach to the construction of the Homotopy equation. Hence, He applied HPM to solve the Lighthill equation [17], Duffing equation [22], and Blasius equation [23], then the idea goes through and has been applied to solve nonlinear wave equations [17], boundary value problems [20]. Babolian et al. [24] applied the homotopy perturbation method to solve the Burgers, the modified Korteweg-de Vries, and regularized long-wave equations.

On the other hand, the HPM has more improved and developed by many engineers and scientists, for instance, a couple of the Laplace transforms and Homotopy perturbation method was implemented by El-Dib et al. [25]. The HPM with two expanding parameters that efficient for some partial nonlinear equations was suggested by He [26] and El-Dib [27]. Also, El-Dib [28] introduced a modified version of the HPM via the multiple scales technique. This new modification works particularly well for the nonlinear oscillators. Furthermore, away from the traditional HPM, Ren et al. [29] made another couple of the HPM and multiple time scales to become a powerful mathematical tool for many nonlinear equations. They displayed that the present procedure may be further afflicted by incorporating several known technologies. It provides solutions to nonlinear equations, whilst the classical perturbation technique became unsuccessful. Moreover, Rabbani [30] introduced a new homotopy perturbation approach for solving main non-linear models through the projection method. A new homotopy perturbation technique for solving linear and nonlinear Schrödinger equations has been addressed by Ayati et al. [31]. Further, by utilizing the HPM, a novel approach in examining the nonlinear Rayleigh-Taylor instability is conducted by El-Dib et al. [32]. Recently, a periodic solution of the cubic nonlinear Klein–Gordon equation using the He-multiple-scales method has been investigated by El-Dib [33]. Also, El-Dib et al. [34] investigated the impact of fractional derivative properties on the periodic analytic solution of the nonlinear oscillations using the HPM. Moreover, He’s-multiple-scale scale to analyze the cubic-quintic Duffing equation has been analyzed by El-Dib et al. [35]. El-Dib [36] presented a stability approach of a fractional-delayed Duffing oscillator. A Nonlinear Instability of a Cylindrical Interface between two MHD Darcian flows has been studied by Moatimid et al. [37]. Further, El-Dib [38] introduced a modified multiple scale technique for the stability of the fractional delayed nonlinear oscillator. Besides, a periodic solution of the fractional sine-Gordon equation has been studied by Shen et al. [39]. Elgazery [40] applied the HPM to give a periodic solution of the Newell-Whitehead-Segel model. Further, for more very useful modification of the homotopy perturbation approach, Yu et al. [41] introduced HPM with an auxiliary parameter for nonlinear oscillators. Also, HPM for Fangzhu oscillator has been used by He et al. [42]. Finally, for more very useful modification of the homotopy perturbation approach, see [43–45].

Motivated by possibility applications in engineering, biology, and physics, which is based on studying the solution of the damped Mathieu equation. Hence, in the present work, our objective is to apply the HPM to linear or nonlinear equations having periodic coefficients such Mathieu equation which has been of great importance among researchers. For the presentation of this article; the rest of the manuscript is systematized as follows: Section 2 is introducing the HPM for the mathematical formulation. The modulation procedure, in detail, is displayed in Section 3. The non-resonance case, stability analysis of the non-resonance case, the resonance case of Ω is near ω , stability analysis of the linear Mathieu equation, and stability analysis for the nonlinear case are represented through Sections 4 to 8, respectively. Finally, the main obtained outcomes are summarized as concluding remarks in Section 9.

2  Mathematical Formulation

To explain the proposed technique, consider the following parametric pendulum equation as an illustrative example:

d2ydt2+μdydt+(ω2+2qcos⁡2Ωt)sin⁡y=0, (1)

where μ,ω,q and Ω are real physical constants. μ is the damping coefficient parameter, ω2 is the square of the natural frequency of the unforced system, q is the amplitude of the parametric excitation, and Ω is the frequency of the parametric excitation. By using the Taylor expansion up to cubic power, the above parametric pendulum equation becomes the following damping cubic nonlinear Mathieu equation:

d2ydt2+μdydt+(ω2+2qcos⁡2Ωt)y=(ω23!+2q3!cos⁡2Ωt)y3+.... (2)

The homotopy perturbation method can be considered as a combination of the classical perturbation technique and the homotopy (whose origin is in the topology [46]), but not restricted to the limitations of traditional perturbation methods. For example, this method does require neither small parameter nor linearization and only requires little iteration to obtain accurate solutions [17] and [18].

We define the two parts of Eq. (2) as L(y) and N(y) , where

L(y)=d2ydt2+ω2y, and N(y)=μdydt+2qycos⁡2Ωt−ky3, (3)

where, k=13!(ω2+2qcos⁡2Ωt) .

Construct the homotopy statement as

H(y,ρ)=L(y)+ρN(y)=0;ρ∈[0,1]. (4)

As in He’s a homotopy perturbation method, it is obvious that when ρ=0 , Eq. (4) becomes the harmonic equation; L(y)=0. Thus,

d2y(t)dt2+ω2y(t)=0. (5)

According to linear differential equations theory, the general solution of Eq. (5) is expressed in terms of two linearly independent solutions, say, eiωt and e−iωt . Thus, the composite solutions may be in the form

y(t)=Aeiωt+A¯e−iωt, (6)

where A and its complex conjugate A¯ are arbitrary constants of integration. Eq. (4) becomes the original nonlinear Mathieu Eq. (2) as ρ=1 . For arbitrary the small parameter ρ , the solution of Eq. (4) can be sought in terms of ρ so that the function y(t) becomes y(t,ρ) . Accordingly, Eq. (4) can rewrite as

H(y,ρ)=(d2dt2+ω2)y(t,ρ)+ρ(μddt+2qcos⁡2Ωt−ky2(t,ρ))y(t,ρ)=0. (7)

It can be noticed that the homotopy function (7) is essentially the same as (4), except for the function y(t,ρ) , which contains embedded the homotopy parameter ρ . The introduction of that parameter within the differential equation is a strategy to redistribute the periodic part between the successive iterations of the homotopy method, and thus increase the probabilities of finding the sought solution. Thus, as ρ moves from 0 to 1, the function y(t,ρ) moves from y0(t) to yapp(t) . Expand the function y(t,ρ) as a power series in the small parameter ρ such that

y(t,ρ)=y0(t)+ρy1(t)+ρ2y2(t)+..., (8)

where yn(t);n=1,2,3,... are unknowns in needs to determine. Substituting this expansion into the homotopy Eq. (7) and equating terms with identical powers of ρ , leads to

ρ0:d2y0dt2+ω2y0=0, (9)

ρ1:d2y1dt2+ω2y1=−μdy0dt−2qcos⁡2Ωty0+ky03, (10)

ρ2:d2y2dt2+ω2y2=−μdy1dt−2qcos⁡2Ωty1+3ky02y1,.... (11)

Eq. (9) can be satisfied by

y0(t)=Aeiωt+A¯e−iωt. (12)

Substituting (12) into Eq. (10) gets the form

d2y1dt2+ω2y1=−iωμ(Aeiωt−A¯e−iωt)−qA(ei(ω+2Ω)t+ei(ω−2Ω)t)−qA¯(e−i(ω+2Ω)t+e−i(ω−2Ω)t)+k(A3e3iωt+3A2A¯eiωt+3A¯2Ae−iωt+A¯3e−3iωt). (13)

Before analyzing the first-order problem, we must distinguish between the two cases. The case of the frequency Ω doesn’t equal the nature frequency ω (which is known as the non-resonance case). The second one is the specific case when Ω approaches ω (which is known as the resonance case).

For arbitrary frequency Ω , there are secular terms that appear in the Eq. (13). Elimination such secular terms require that the arbitrary constant A be zero. This means that the expansion (8) cannot be successful to obtain a valid solution for exciting homotopy Eq. (7).

3  The Modulation Procedure

To obtain uniform expansions for problems of this kind, the expansion (8) needs to be modified. If we modulate the initial solution (6) so that the constant A becomes A(τ) with τ=ρt, such that

dAdt=ρdAdτ and d2Adt2=ρ2dAdτ2. (14)

Then Eq. (12) in the modulate case becomes

Y0(t,τ)=A(τ)u0(t)+A¯(τ)u¯0(t), (15)

where

u0(t)=eiωt and u¯0(t)=e−iωt. (16)

Consequently, the homotopy state, Eq. (7), in the modulated form becomes

(d2dt2+ω2)Y(t,τ,ρ)+ρ(μddt+2qcos⁡2Ωt)Y(t,τ,ρ)=ρkY3(t,τ,ρ). (17)

It is convenient to choose the modulated function Y(t,τ,ρ) in separated variables as

Y(t,τ,ρ)=A(τ)u(t,ρ)+A¯(τ)u¯(t,ρ). (18)

The function u(t,ρ) can be expanded as a power series in the small parameter ρ such that

u(t,ρ)=u0(t)+ρu1(t)+ρ2u2(t)+..., (19)

where un(t);n=1,2,3,... are unknowns to be evaluated. If the expansion (19) is substituted into Eq. (18) then gets

Y(t,τ,ρ)=A(τ)(u0(t)+ρu1(t)+ρ2u2(t)+...)+cc.=Y0(t,τ)+ρY1(t,τ)+ρ2Y2(t,τ)+..., (20)

where cc. indicates to the complex conjugate for the preceding terms and

Yn(t,τ)=A(τ)un(t)+A¯(τ)u¯n(t). (21)

It is noted that:

ddtY(t,τ,ρ)=ddt[A(τ)u(t,ρ)+cc]=A(τ)u˙(t,ρ)+ρu(t,ρ)A′(τ)+cc, (22)

and

d2dt2Y(t,τ,ρ)=A(τ)u¨(t,ρ)+2ρA′(τ)u˙(t,ρ)+ρ2u(t,ρ)A′′(τ)+cc, (23)

where dots indicate differentiation concerning the time t , while dashes refer to the derivative for the time modulate τ. Substituting (18) into Eq. (17) using (22) and (23) gives

A(u¨+ω2u)+ρ(2u˙A′+μAu˙+2qAucos⁡2Ωt)+ρ2(A′′+μA′)u−ρk(A3u3+3A2A¯u2u¯)+cc.=0. (24)

Eq. (24) remains to obey the same homotopy concept because it’s become the same harmonic Eq. (5) as ρ→0 . limρ→1⁡A′=ddτ(limρ→1⁡A)=0, besides, consequently the original Eq. (2) is found.

In the light of Eq. (19), the modulate homotopy Eq. (24) will be expanded as a power series in ρ so that the following non-homogenous harmonic equations are imposed

ρ0:A(u¨0+ω2u0)+A¯(u¨¯0+ω2u¯0)=0, (25)

ρ1:A(u¨1+ω2u1)+2u˙0A′+μAu˙0+2qAu0cos⁡2Ωt−kA2(Au03+3A¯u02u¯0)+cc.=0, (26)

ρ2:A(u¨2+ω2u2)+2u˙1A′+μAu˙1+2qAu1cos⁡2Ωt+u0A′′+μA′u0−3kA2[Au02u1+A¯(u02u¯1+2u0u¯0u1)]+cc.=0. (27)

It is noted that Eq. (25) has been satisfied by Eq. (16) and the zero-order solution for Eq. (17) as approved in Eq. (15). Substituting Eq. (16) into Eq. (26) becomes

A(u¨1+ω2u1)+[iω(2A′+μA)−3kA2A¯]eiωt+qA(ei(ω+2Ω)t+ei(ω−2Ω)t)−kA3e3iωt+cc.=0. (28)

This equation contains secular terms at the non-resonance case and other secular terms when the applied frequency Ω approaches the natural frequency ω .

4  The Non-Resonance Case

The analysis in this case concerned with the arbitrary chosen for the applied frequency Ω , in Eq. (28). At this stage, secular terms are removed when

A′+12μA+3ik2ωA2A¯=0, (29)

with its complex conjugate one. This leads to obtaining the valid function u1(t) as

u1(t)=q4Ω(ei(ω+2Ω)t(Ω+ω)+ei(ω−2Ω)t(Ω−ω))+kA28ω2e3iωt. (30)

Consequently, the solution of the first-order problem is formulated as

Y1(t,τ)=qA(τ)4Ω(ei(ω+2Ω)t(ω+Ω)+ei(ω−2Ω)t(ω−Ω))+qA¯(τ)4Ω(e−i(ω+2Ω)t(ω+Ω)+e−i(ω−2Ω)t(ω−Ω))+k8ω2(A3(τ)e3iωt+A¯3(τ)e−3iωt). (31)

Substituting Eqs. (16) and (30) into Eq. (27), using Eq. (29), yields

A(u¨2+ω2u2)+[A′′+μA′+q2A2(Ω2−ω2)−3k28ω2A3A¯2]eiωt+3k2A4A¯8ω2e3iωt−3k2A58ω2e5iωt+q2A4Ω(Ω2−ω2)[(Ω−ω)ei(ω+4Ω)t+(Ω+ω)ei(ω−4Ω)t]+3kqA2A¯2ω(Ω2−ω2)[(Ω−2ω)ei(ω+2Ω)t−(Ω+2ω)ei(ω−2Ω)t]−3kqA34Ω(Ω2−ω2)[(Ω−ω)ei(3ω+2Ω)t+(Ω+ω)ei(3ω−2Ω)t]+cc.=0. (32)

The valid solution requires to be removed the terms that producing unbounded solution. These terms imply the following nonlinear solvability condition:

A′′+μA′+q2A2(Ω2−ω2)−3k28ω2A3A¯2=0. (33)

The second-order solution is found to be

Y2(t,τ)=q2A(τ)32Ω2(Ω2−ω2)[(Ω−ω)(ω+2Ω)ei(ω+4Ω)t+(Ω+ω)(ω−2Ω)ei(ω−4Ω)t]−3kqA3(τ)16Ω(Ω2−ω2)2[(Ω−ω)2(2ω+Ω)ei(3ω+2Ω)t+(Ω+ω)2(Ω−2ω)ei(3ω−2Ω)t]+3kqA2(τ)A¯(τ)8Ω2ω(Ω2−ω2)[(Ω−2ω)(ω+Ω)ei(ω+2Ω)t+(Ω+2ω)(Ω−ω)ei(ω−2Ω)t]−k2A5(τ)64ω2e5iωt+3k2A4(τ)A¯(τ)64ω2e3iωt+cc. (34)

If the accuracy to the second-order perturbation is enough, then the approximate solution at the non-resonance case is formulated by substituting Eqs. (15), (16), (31) and (34) into Eq. (20), and setting ρ=1 , gets

Y(t)=limρ→1,τ→t⁡Y(t,τ,ρ)=A(t)eiωt+k8ω2A3(t)e3iωt−k2A5(t)64ω2e5iωt+3k2A4(t)A¯(t)64ω2e3iωt+qA(t)4Ω(ei(ω+2Ω)t(ω+Ω)+ei(ω−2Ω)t(ω−Ω))+3kqA2(t)A¯(t)8Ω2ω(Ω2−ω2)[(Ω−2ω)(ω+Ω)ei(ω+2Ω)t+(Ω+2ω)(Ω−ω)ei(ω−2Ω)t]+q2A(t)32Ω2(Ω2−ω2)[(Ω−ω)(ω+2Ω)ei(ω+4Ω)t+(Ω+ω)(ω−2Ω)ei(ω−4Ω)t]−3kqA3(t)16Ω(Ω2−ω2)2[(Ω−ω)2(2ω+Ω)ei(3ω+2Ω)t+(Ω+ω)2(Ω−2ω)ei(3ω−2Ω)t]+cc. (35)

5  Stability Analysis for the Non-Resonance Case

The stability criteria in the non-resonance case can be obtained from solving Eq. (29). One may use the following polar form [16]:

A(τ)=12ξ(τ)eiη(τ), (36)

with real the unknown functions ξ(τ) and η(τ) . Insert Eq. (36) into the first-order solvability condition (29) which will separate into real and imaginary parts and gives

ξ(τ)=ξ0e−12μτ and η(τ)=3k2μωξ0e−12μτ+η0, (37)

where, ξ0 and η0 are integration constants. The stability criteria in the non-resonance case require that μ>0.

6  The Resonance Case Ω is Near ω

Return to the first-order problem Eq. (28) and re-analyzed it because of the nearness of Ω to ω . We express this approach by introducing the detuning parameter σ [16] such that

Ω=ω+ρσ. (38)

Accordingly, we have

−i(ω−2Ω)t=iωt+2iστ. (39)

Elimination of secular terms from Eq. (28), because of Eqs. (38) and (39) yields

A′+12μA−iq2ωA¯e2iστ+3ik2ωA2A¯=0. (40)

The first-order solution in this case is

Y1(t,τ)=q4Ω(ω+Ω)(Aei(ω+2Ω)t+A¯e−i(ω+2Ω)t)+k8ω2(A3e3iωt+A¯3e−3iωt). (41)

Using Eq. (41) with Eq. (27), we obtain the uniform solution for the second-order problem, and the following solvability is presented:

A′′+μA′+q24Ω(ω+Ω)A+kq8ω2A3e−2iστ−3kq4Ω(ω+Ω)AA¯2e2iσt−3k28ω2A3A¯2=0, (42)

with its complex conjugate. The valid function Y2(t,τ) is given by

Y2(t,τ)=3ik64ω3(2A′+μA+2ikA2A¯)A2e3iωt−3k2192ω4A5e5iωt+q232Ω2(ω+Ω)(ω+2Ω)Aei(ω+4Ω)t+iq8Ω2(ω+Ω)2[(A′+12μA)(ω+2Ω)+3ikA2A¯]ei(ω+2Ω)t−kq16(2ω+Ω)(ω+Ω)[3Ω(ω+Ω)−12ω2]A3ei(3ω+2Ω)t+cc. (43)

The approximate solution up to the second-order is formulated by substituting from Eqs. (15), (16), (41) and (43) into Eq. (16) gets

Y(t)=limρ→1τ→t⁡(Y0+ρY1+ρ2Y2)=Aeiωt+A¯e−iωt+k8ω2[A+3i8ω(2A′+μA+2ikA2A¯)]A2e3iωt−3k2192ω4A5e5iωt+q4Ω(ω+Ω){A+iq2Ω(ω+Ω)[(A′+12μA)(ω+2Ω)+3ikA2A¯]}ei(ω+2Ω)t+q232Ω2(ω+Ω)(ω+2Ω)Aei(ω+4Ω)t−kq16(2ω+Ω)(ω+Ω)[3Ω(ω+Ω)−12ω2]A3ei(3ω+2Ω)t+cc. (44)

7  Stability Analysis of the Linear Mathieu Equation

In the limiting case as k→0 into Eq. (2), linear damping Mathieu equation arrived. In this case, the two solvability conditions (40) and (42) that produced at the resonance case of Ω is near ω having the following limit case:

A′+12μA−iq2ωA¯e2iστ=0, (45)

A′′+μA′+q24Ω(ω+Ω)A=0. (46)

The first-order solvability condition (45) can be used to find the stability picture at the resonance case. The second-order solvability condition (46) can be used to find the value of the detuning parameter σ .

It is easy to show that the Eq. (45) can be satisfied by the form

A(τ)=[(σ+q2ω)sin⁡Θτ+iΘcos⁡Θτ]e(iσ−12μ)τ, (47)

where, the parameter μ must be positive, to find a damping solution. The argument Θ is given by the following characteristic equation:

Θ2=σ2−q24ω2. (48)

The parameter σ can be evaluated by substituting Eq. (47) into the second-order solvability condition (46) to gets

σ=−q2ωorσ=±ω4q(μ2−q2ω2−q2Ω(Ω+ω)). (49)

The use of the first value of σ , Eq. (48) yields a zero solution for Eq. (45). For a non-zero solution, the other values σ are conforming. Inserting Eq. (49) into Eq. (48) gets

Θ2=ω216q2(μ2−q2ω2−q2Ω(Ω+ω))2−q24ω2. (50)

The stability criteria require that the right-hand-side of Eq. (50) be positive, which implies that

(μ2−q2ω2−q2Ω(Ω+ω))2−4q4ω4>0. (51)

Stability condition (51) can be rearranged in powers of the applied frequency Ω as

Ω2(μ2ω2−3q2)+Ωω(μ2ω2−3q2)−q2ω2>0, (52)

and

Ω2(μ2ω2+q2)+Ωω(μ2ω2+q2)−q2ω2<0. (53)

The transition curves separating stable state from an unstable state corresponding to

Ω1=−ω(μ2ω2−3q2)−ωμ4ω4−2μ2ω2q2−3q42(μ2ω2−3q2), (54)

and

Ω2=−ω(μ2ω2+q2)−ω(μ2ω2+q2)(μ2ω2+5q2)2(μ2ω2+q2). (55)

8  Stability Analysis for the Nonlinear Case

The first-order solvability condition (40) can be used to find the stability picture at the resonance case. The second-order solvability condition (42) can be used to find the value of the detuning parameter σ .

To relax the periodic term into the Eq. (40) we let

A(τ)=[α(τ)+iβ(τ)]eiστ, (56)

with real functions α and β . Insert Eq. (56) into Eq. (40), separating real and imaginary parts yields:

α′+12μα−(σ+q2ω+3k2ω(α2+β2))β=0, (57)

β′+12μβ+(σ−q2ω−3k2ω(α2+β2))α=0. (58)

In order to solve the above coupled nonlinear Eqs. (57) and (58), we may discuss the behavior at the steady-state response. This case is corresponding to the case of d..dτ=0 . If the solutions of Eqs. (57) and (58), at the steady-state, are represented by α0 and β0 , which are given by

12μα0−(σ+q2ω+3kr22ω)β0=0, (59)

12μβ0+(σ−q2ω−3kr22ω)α0=0, (60)

where r2=α02+β02 is used. Eqs. (59) and (60) are two coupled algebraic equations in α0 and β0 . For nontrivial solutions in α0 and β0 , we obtain

σ2=(q+3kr2)4ω22−14μ2. (61)

Besides, the constants α0 and β0 may be chosen as

α0=(σ+q+3kr22ω), and β0=12μ. (62)

Squaring both equations in (62) and adding we get

(σ+q+3kr22ω)2=r2−14μ2. (63)

Combing Eq. (61) with Eq. (63) yields

σ=2ω2r2−(q+3kr2)22ω(q+3kr2). (64)

In order to find a constrain for a bounded solution we may modulate the functions α and β as

α(τ)=α0+α1(τ)andβ(τ)=β0+β1(τ), (65)

where the functions α1(τ) and β1(τ) refer to a small deviation from the steady-state solution α0 and β0 . Then the system of Eqs. (57) and (58) in the linearizing form becomes

α1′+(12μ−3kωα0β0)α1−(σ+q+3kr22ω+3kωβ02)β1=0, (66)

β′+(12μ−3kωα0β0)β1+(σ−q+3kr22ω−3kωα02)α1=0. (67)

The above system is two coupled linear differential equations of first-order in the two functions α1 and β1 . This system can be satisfied by

α1(τ)=(σ+q+3kr22ω+3kωβ02)e−(12μ−3kωα0β0)τsin⁡Θτ, (68)

β1(τ)=e−(12μ−3kωα0β0)τΘcos⁡Θτ, (69)

where Θ is given by the following characteristic equation:

Θ2=(σ−q+9kr22ω+3k4ωμ2)(σ+q+3kr22ω+3k4ωμ2). (70)

where relations (62) are used. This characteristic equation depends on the two related parameters σ and r2 . This relation between them is given in Eq. (61) or in Eq. (64).

With the help of the second-order solvability condition (42) one can find an expression for both the unknowns σ and r2 in terms of the frequency Ω . To accomplish this, one may substitute the steady-state solution

A(τ)=(α0+iβ0)eiστ, (71)

into the second-order solvability condition (42). Separating the real and imaginary parts, produces the following relations, between the parameters σ , Ω and r2 :

σ2−q24Ω(ω+Ω)+116kqμ2(Ω2+ωΩ−6ω2ω2Ω(ω+Ω))−[kq8(Ω2+ωΩ−6ω2ω2Ω(ω+Ω))−3k28ω2(σ+q+3kr22ω)]r2=0, (72)

σ+kq8(Ω2+ωΩ+6ω2ω2Ω(ω+Ω))(σ+q+3kr22ω)+3k22ω2r2=0, (73)