Construction of Approximate Analytical Solutions to Strongly Nonlinear Coupled van der Pol Oscillators

Abstract

Using nonlinear theory to research vibration model of engineering system has important theoretical and practical significance. Multi-degree-of-freedom (MDOF) coupled van der Pol oscillator is a typical model in the nonlinear vibration; many complex dynamic problems in practical engineering can be simplified as this model to be solved in the end. This paper discusses a class of two-degrees-of-freedom (2-DOF) coupled van der Pol oscillator, which was divided into three parameters of different situations α₁≠α₂, β₁≠β₂, and γ₁≠γ₂ to discuss. Employing symbolic software such as Mathematica for those problems, the explicit analytical solutions of frequency ω and displacements x₁(t) and x₂(t) are well formulated. Results showed that the homotopy analysis method (HAM) can effectively deal with this kind of parameter of different coupled vibrators, just request the values of some parameters are not too big. Finally, we got four important theorems to simplify the solution of the nonlinear system.

1. Introduction

Most of practical engineering problems can be governed by multi-degree-of-freedom (MDOF) nonlinear systems. The significance of MDOF nonlinear dynamical systems is mainly due to its complex dynamics, global bifurcations, and chaotic motions; the intensive research subjects are thus at the forefront of nonlinear dynamics. Due to lack of analytical tools and methods to study MDOF nonlinear systems, it is extremely challenging to develop the theories of MDOF nonlinear systems and to give systematic applications to engineering problems.

Enlightened from the basic concepts of the homotopy in topology, Liao [1, 2] proposed an analytical method for solving a kind of nonlinear problems, namely, the homotopy analysis method (HAM). Unlike the classical perturbation techniques, the clear-cut feature of the HAM is independent of small parameters in the governing equations of motion. Inasmuch as the HAM can overcome the foregoing restrictions of conventional perturbation techniques, this method has been widely applied to solve a variety of nonlinear problems [3 –7].

The studies of the limit cycles (LC) of the van der Pol oscillator have attracted a considerable number of authors [8, 9]. Odani obtained the LC of the van der Pol equation that is not algebraic in [10]. Wirkus and Rand [11] discussed the dynamics of two coupled van der Pol oscillators with delay coupling. Hirano and Rybicki [12] presented the existence of limit cycles for coupled van der Pol equations. Bi [13] studied the dynamical analysis of two coupled parametrically excited van der Pol oscillators. Rompala et al. [14] proposed the dynamics of three coupled van der Pol oscillators; they investigated the existence of the in-phase mode in which the first two oscillators have the same behavior. Recently, Chen and Liu [15, 16] studied the LC of van der Pol and van der Pol-Duffing equation for single degree of freedom (SDOF) by the HAM. Kuznetsov et al. [17] obtained the phase dynamics and structure of synchronization tongues of coupled van der Pol-Duffing oscillators. Kimiaeifar et al. [18] made use of HAM to give the analytical solution for van der Pol-Duffing oscillators for SDOF. Qian et al. [19] derived the highly convergent solutions for MDOF coupled van der Pol oscillators. Liao [20] elucidated some new definitions and theorems of the HAM. As a result of the existence of nonlinear terms, we must study different kinds of systems, no matter the physical significance and complication of calculation. With HAM as an analytical approximation method, we can verify its effectiveness through specific cases. In this paper we study a pair of van der Pol oscillators in which there are cubic couple terms. Because of its flexibility, the present techniques can also be further generalized to analyze more complicated nonlinear MDOF dynamical systems that can only be analyzed by numerical methods.

In this paper, the HAM is applied to study the accurate approximate analytical solutions for a class of coupled van der Pol oscillators. We discuss a 2-DOF coupled van der Pol system in three circumstances with different parameters α₁≠α₂, β₁≠β₂, and γ₁≠γ₂. The validity and accuracy of the HAM are also performed. The auxiliary parameter ℏ is firstly modified to ensure the convergence of solution series. Also, the relationship of the frequency, amplitudes, and parameters is investigated through figures. Then we acquire the approximate value of ω and expressions of x₁(t) and x₂(t). The comparison between approximate solutions and exact ones is conducted simultaneously. It is shown that the periodic solutions of the HAM are in excellent agreement with the numerical integration ones, even if time t progresses to a certain large domain. In what follows, Section 2 presents the HAM of the 2-DOF coupled van der Pol system. In Section 3, numerical comparisons are carried out to authenticate the correctness and accuracy of the present method. Moreover, we got four important theorems to simplify the solution of the nonlinear system in Section 4. Finally, the paper ends with concluding remarks in Section 5.

2. Homotopy Analysis Method

In this section, we consider the following types of the 2-DOF strongly nonlinear coupled van der Pol oscillators

{\ddot{x}}_{1} + α_{1} (x_{1}^{2} - 1) {\dot{x}}_{1} + β_{1} x_{1} + δ_{1} x_{1}^{3} + γ_{1} (x_{1}^{3} - x_{2}^{3}) = 0,

(1a)

{\ddot{x}}_{2} + α_{2} (x_{2}^{2} - 1) {\dot{x}}_{2} + β_{2} x_{2} + δ_{2} x_{2}^{3} + γ_{2} (x_{2}^{3} - x_{1}^{3}) = 0,

(1b)

where x₁(t) and x₂(t) are unknown functions, α_i, β_i, γ_i, and δ_i (i = 1, 2) are the known physical parameters, and x₁ and x₂ are mutual coupling.

Introduce a new variable τ = ωt and substitute it into the systems (1a) and (1b); we get

ω^{2} x_{1}^{′′} + ω α_{1} (x_{1}^{2} - 1) x_{1}^{'} + β_{1} x_{1} + δ_{1} x_{1}^{3} + γ_{1} (x_{1}^{3} - x_{2}^{3}) = 0,

(2a)

ω^{2} x_{2}^{′′} + ω α_{2} (x_{2}^{2} - 1) x_{2}^{'} + β_{2} x_{2} + δ_{2} x_{2}^{3} + γ_{2} (x_{2}^{3} - x_{1}^{3}) = 0,

(2b)

with the initial conditions

x_{1} (0) = a, x_{1}^{'} (0) = b, x_{2} (0) = c, x_{2}^{'} (0) = 0 .

(3)

Apostrophe is on behalf of the new variable's differential.

Assuming the periodic solution of systems (2a) and (2b) can be represented by the basis function {cos⁡(kτ), sin(kτ) ∣ k = 0, 1, 2, …}

x_{1} (τ) = \sum_{k = 0}^{+ \infty} (a_{1, k} \cos k τ + b_{1, k} \sin k τ),

(4a)

x_{2} (τ) = \sum_{k = 0}^{+ \infty} (a_{2, k} \cos k τ + b_{2, k} \sin k τ) .

(4b)

For the initial guess solution, we set

x_{1, 0} (τ) = a_{0} \cos τ + b_{0} \sin τ, x_{2, 0} (τ) = c_{0} \cos τ,

(5)

and the auxiliary linear operator is expressed as

L (\begin{pmatrix} Φ_{1} (τ, q) \\ Φ_{2} (τ, q) \end{pmatrix}) = ω_{0}^{2} (\begin{pmatrix} \frac{\partial^{2} Φ_{1} (τ, q)}{\partial τ^{2}} + Φ_{1} (τ, q) \\ \frac{\partial^{2} Φ_{2} (τ, q)}{\partial τ^{2}} + Φ_{2} (τ, q) \end{pmatrix}) .

(6)

According to the systems (2a) and (2b), nonlinear operator is defined as follows:

N (\begin{pmatrix} Φ_{1} (τ, q) \\ Φ_{2} (τ, q) \end{pmatrix}) = (\begin{matrix} Ω^{2} (q) \frac{\partial^{2} Φ_{1} (τ, q)}{\partial τ^{2}} + Ω (q) α_{1} [Φ_{1}^{2} (τ, q) - 1] \frac{\partial Φ_{1} (τ, q)}{\partial τ} + β_{1} Φ_{1} (τ, q) \\ + δ_{1} Φ_{1}^{3} (τ, q) + γ_{1} [Φ_{1}^{3} (τ, q) - Φ_{2}^{3} (τ, q)] \\ Ω^{2} (q) \frac{\partial^{2} Φ_{2} (τ, q)}{\partial τ^{2}} + Ω (q) α_{2} [Φ_{2}^{2} (τ, q) - 1] \frac{\partial Φ_{2} (τ, q)}{\partial τ} + β_{2} Φ_{2} (τ, q) \\ + δ_{2} Φ_{2}^{3} (τ, q) + γ_{2} [Φ_{2}^{3} (τ, q) - Φ_{1}^{3} (τ, q)] \end{matrix}) .

(7)

Taking H(t) = 1, the zero-order deformation equation can be written as

(1 - q) L (\begin{pmatrix} Φ_{1} (τ, q) - x_{1, 0} (τ) \\ Φ_{2} (τ, q) - x_{2, 0} (τ) \end{pmatrix}) = q ℏ N (\begin{pmatrix} Φ_{1} (τ, q) \\ Φ_{2} (τ, q) \end{pmatrix}),

(8)

with the initial conditions

Φ_{1} (0, q) = a (q), {\frac{\partial Φ_{1} (τ, q)}{\partial τ} |}_{τ = 0} = b (q),

(9a)

Φ_{2} (0, q) = c (q), {\frac{\partial Φ_{2} (τ, q)}{\partial τ} |}_{τ = 0} = 0 .

(9b)

When q = 0, (8), (9a), and (9b) have solutions

\begin{matrix} Φ_{1} (τ, 0) = x_{1, 0} (τ), Φ_{2} (τ, 0) = x_{2, 0} (τ), \\ Ω (0) = ω_{0} . \end{matrix}

(10)

When q = 1, the zero-order deformation equations (8), (9a), and (9b) are equivalent to the original equations (2a), (2b), and (3), so

\begin{matrix} Φ_{1} (τ, 1) = x_{1} (τ), Φ_{2} (τ, 1) = x_{2} (τ), \\ Ω (1) = ω . \end{matrix}

(11)

When the embedded variable q increases from 0 to 1, Φ_i (τ,q) changes from the initial guess solution x_i,0(τ) to an unknown exact solution x_i(τ) (i = 1, 2). Likewise, Ω(q) changes from an initial guess frequency ω₀ to the physical frequency ω.

Using Taylor expansion theorem and (10), we get

Φ_{1} (τ, q) = x_{1, 0} (τ) + \sum_{m = 1}^{+ \infty} x_{1, m} (τ) q^{m},

(12a)

Φ_{2} (τ, q) = x_{2, 0} (τ) + \sum_{m = 1}^{+ \infty} x_{2, m} (τ) q^{m},

(12b)

Ω (q) = ω_{0} + \sum_{m = 1}^{+ \infty} ω_{m} q^{m},

(12c)

{x_{1, m} (τ) = \frac{1}{m!} \frac{\partial^{m} Φ_{1} (τ, q)}{\partial q^{m}} |}_{q = 0},

(13a)

{x_{2, m} (τ) = \frac{1}{m!} \frac{\partial^{m} Φ_{2} (τ, q)}{\partial q^{m}} |}_{q = 0},

(13b)

{ω_{m} = \frac{1}{m!} \frac{\partial^{m} Ω (q)}{\partial q^{m}} |}_{q = 0} .

(13c)

Assuming that the auxiliary parameter ℏ is selected appropriately, (12a), (12b), and (12c)'s power series are converging when q = 1. By (13a), (13b), and (13c), we get the series solution

x_{1} (τ) = x_{1, 0} (τ) + \sum_{m = 1}^{+ \infty} x_{1, m} (τ),

(14a)

x_{2} (τ) = x_{2, 0} (τ) + \sum_{m = 1}^{+ \infty} x_{2, m} (τ),

(14b)

ω = ω_{0} + \sum_{m = 1}^{+ \infty} ω_{m} .

(14c)

For convenience, we define the vector as follows:

x_{1, n} = {x_{1, 0} (τ), x_{1, 1} (τ), \dots, x_{1, n} (τ)},

(15a)

x_{2, n} = {x_{2, 0} (τ), x_{2, 1} (τ), \dots, x_{2, n} (τ)},

(15b)

ω_{n} = {ω_{0}, ω_{1}, \dots, ω_{n}} .

(15c)

By differentiating the zeroth-order deformation equation (8) m times with respect to q, and both sides of the resulting equation are divided by m!, and assume q = 0; we obtain m order deformation equation

L (\begin{pmatrix} x_{1, m} (τ) - χ_{m} x_{1, m - 1} (τ) \\ x_{2, m} (τ) - χ_{m} x_{2, m - 1} (τ) \end{pmatrix}) = ℏ (\begin{pmatrix} R_{1, m} (x_{1, m - 1}, ω_{m - 1}) \\ R_{2, m} (x_{2, m - 1}, ω_{m - 1}) \end{pmatrix}),

(16)

with the initial conditions

\begin{matrix} x_{1, m} (0) = a_{m}, x_{1, m}^{'} (0) = b_{m}, x_{2, m} (0) = c_{m}, x_{2, m}^{'} (0) = 0, (m \geq 1), \\ (\begin{pmatrix} R_{1, m} (x_{1, m - 1}, ω_{m - 1}) \\ R_{2, m} (x_{2, m - 1}, ω_{m - 1}) \end{pmatrix}) = \frac{1}{(m - 1)!} {\frac{\partial^{m - 1} N (\begin{smallmatrix} Φ_{1} (τ, q) \\ Φ_{2} (τ, q) \end{smallmatrix})}{\partial q^{m - 1}} |}_{q = 0} . \end{matrix}

(17)

According to the solution expression principles and the definition of linear operator L, by (16), we obtained solutions x_1,m(τ) and x_2,m(τ) without secular terms τsinτ and τcos⁡τ. The right-hand side of (16) cannot contain sinτ, cos⁡τ, or similar items. We set the coefficient of these items are 0, and we get

\frac{1}{π} \int_{0}^{2 π} R_{1, m} (x_{1, m - 1}, ω_{m - 1}) \cos τ d τ = 0,

(18a)

\frac{1}{π} \int_{0}^{2 π} R_{2, m} (x_{2, m - 1}, ω_{m - 1}) \cos τ d τ = 0,

(18b)

\frac{1}{π} \int_{0}^{2 π} R_{1, m} (x_{1, m - 1}, ω_{m - 1}) \sin τ d τ = 0,

(18c)

\frac{1}{π} \int_{0}^{2 π} R_{2, m} (x_{2, m - 1}, ω_{m - 1}) \sin τ d τ = 0 .

(18d)

Then

R_{1, m} (x_{1, m - 1}, ω_{m - 1}) = \sum_{n = 1}^{φ_{1} (m)} [a_{m, n}^{1} \cos n τ + b_{m, n}^{1} \sin n τ],

(19a)

R_{2, m} (x_{2, m - 1}, ω_{m - 1}) = \sum_{n = 1}^{φ_{2} (m)} [a_{m, n}^{2} \cos n τ + b_{m, n}^{2} \sin n τ],

(19b)

a_{m, n}^{1} = \frac{1}{π} \int_{0}^{2 π} R_{1, m} (x_{1, m - 1}, ω_{m - 1}) \cos (n τ) d τ,

(20a)

a_{m, n}^{2} = \frac{1}{π} \int_{0}^{2 π} R_{2, m} (x_{2, m - 1}, ω_{m - 1}) \cos (n τ) d τ,

(20b)

b_{m, n}^{1} = \frac{1}{π} \int_{0}^{2 π} R_{1, m} (x_{1, m - 1}, ω_{m - 1}) \sin (n τ) d τ,

(20c)

b_{m, n}^{2} = \frac{1}{π} \int_{0}^{2 π} R_{2, m} (x_{2, m - 1}, ω_{m - 1}) \sin (n τ) d τ .

(20d)

Using (16), (18a), (18b), (18c), and (18d), we can gradually get ω_k, a_k, b_k, and c_k (k = 0, 1, 2, …).

3. Numerical Simulation

In this section, we get the periodic solution of systems (1a) and (1b) by HAM under the condition of the given parameters α_i, β_i, γ_i, and δ_i (i = 1, 2). First we consider the case α₁≠α₂; take α₁ = 1, α₂ = 99/100, β₁ = β₂ = 1, γ₁ = γ₂ = 1, δ₁ = δ₂ = 1, and then put them into (18a), (18b), (18c), and (18d); we can get the following nonlinear algebraic equations:

6 a_{0}^{3} - 3 c_{0}^{3} + a_{0}^{2} b_{0} ω_{0} + b_{0} (- 4 + b_{0}^{2}) ω_{0} + a_{0} (4 + 6 b_{0}^{2} - 4 ω_{0}^{2}) = 0,

(21a)

6 b_{0}^{3} - a_{0} (- 4 + a_{0}^{2}) ω_{0} - a_{0} {b_{0}^{2} ω}_{0} + b_{0} (4 + 6 a_{0}^{2} - 4 ω_{0}^{2}) = 0,

(21b)

- 3 a_{0}^{3} - 3 a_{0} b_{0}^{2} + 4 c_{0} + 6 c_{0}^{3} - 4 c_{0} ω_{0}^{2} = 0,

(21c)

100 a_{0}^{2} b_{0} + 100 b_{0}^{3} + 33 c_{0} (- 4 + c_{0}^{2}) ω_{0} = 0 .

(21d)

Solve the equations; we have a₀ = 2, b₀ = 0, c₀ = 2, and ω₀ = 2.

Next make use of the relationship between frequency ω and parameter ℏ to determine the value range of ℏ. Figure 1 is the curve of ω-ℏ by the fourth-order approximation; when − 0.28 < ℏ < − 0.18, the series of ω is convergent. The selection of different auxiliary parameters ℏ can improve the accuracy of the final series solution; for simplicity, take ℏ = − 0.2.

Figure 1:

The fourth-order approximation curve of ω versus ℏ.

We further explore the relationship between the frequency, amplitude, and the parameters. According to (3), the amplitudes of the two components x₁(t), x₂(t) are a, c in periodic solution of the system. In this example, since α₁≠α₂ and other parameter values are the same, the focus of the research is the relationship between ω, a, c and α₁, α₂. First, set α₁ as an unknown parameter, and set α₂ as a known quantity (α₂ = 99/100). According to the HAM, the second-order approximate expression of ω(α₁), a(α₁), and c(α₁) can be written, respectively, as follows:

ω (α_{1}) = \frac{44779}{256000000} α_{1}^{4} - \frac{1731873}{6400000000} α_{1}^{3} - \frac{23516728323}{1280000000000} α_{1}^{2} + \frac{986875659}{102400000000} α_{1} + \frac{78481134113267}{40960000000000},

(22a)

a (α_{1}) = - \frac{121}{4000000} α_{1}^{4} + \frac{714791}{800000000} α_{1}^{3} - \frac{2575699901}{240000000000} α_{1}^{2} + \frac{81116409}{2560000000} α_{1} + \frac{2458329452757}{1280000000000},

(22b)

c (α_{1}) = \frac{13}{115200} α_{1}^{4} - \frac{363}{1280000} α_{1}^{3} - \frac{3893383}{384000000} α_{1}^{2} + \frac{138557859}{12800000000} α_{1} + \frac{6213267}{3200000} .

(22c)

Then set α₂ as an unknown parameter, α₁ as a known quantity (α₁ = 1); the second-order approximate expression of ω(α₂), a(α₂), and c(α₂) can be written, respectively, as follows:

ω (α_{2}) = \frac{1937}{11059200} α_{2}^{4} - \frac{757}{2764800} α_{2}^{3} - \frac{101519}{5529600} α_{2}^{2} + \frac{2989}{307200} α_{2} + \frac{784659}{409600},

(23a)

a (α_{2}) = \frac{17}{207360} α_{2}^{4} + \frac{919}{1036800} α_{2}^{3} - \frac{23207}{1036800} α_{2}^{2} + \frac{3691}{115200} α_{2} + \frac{1159}{600},

(23b)

c (α_{2}) = - \frac{31}{115200} α_{2}^{3} + \frac{19}{12800} α_{2}^{2} + \frac{419}{38400} α_{2} + \frac{44471}{23040} .

(23c)

Figure 2 shows the image of ω(α₁) and ω(α₂) drawn in the same Cartesian coordinate system, namely, the curves of ω-α₁ and ω-α₂.

Figure 2:

The second-order approximation curves of ω-α₁ and ω-α₂.

One can see from Figure 2 that if α₁, α₂ are within the suitable range (− 1/2 < α₁, α₂ < 3/2), and the gap between α₁ and α₂ is very small (|α₁ − α₂| < 1/50), then the approximate value of frequency ω can be stabilized in the vicinity of 1.9. Thus it can be seen that the setting of α₁, α₂ in this example is reasonable.

Figures 3 and 4 give the fourth-order analytical approximate phase diagram and time history curve, respectively. The fourth-order analytical approximation agrees well with the solution obtained by the numerical integration method. Here the initial conditions that numerical integration method meets are x₁(0) = 2, x₁′(0) = 0, x₂(0) = 2, and x₂′(0) = 0.

Figure 3:

Comparison of the phase portrait curves of the fourth-order approximation with the numerical integration method.

Figure 4:

Comparison of the time history responses s of the fourth-order approximation with the numerical integration method.

The specific expressions of the fourth-order analytical approximate that the systems (1a) and (1b) satisfy the above initial conditions are

x_{1} (t) = 1.84957290 \cos ω t + 0.32613610 \sin ω t + 0.09496963 \cos 3 ω t - 0.08934836 \sin 3 ω t - 0.00426108 \cos 5 ω t - 0.01258828 \sin 5 ω t - 0.00140884 \cos 7 ω t + 0.00008743 \sin 7 ω t + 0.00002442 \cos 9 ω t + 0.00009499 \sin 9 ω t,

(24a)

x_{2} (t) = 1.84958869 \cos ω t + 0.32448489 \sin ω t + 0.09453651 \cos 3 ω t - 0.08812457 \sin 3 ω t - 0.00411388 \cos 5 ω t - 0.01230200 \sin 5 ω t - 0.00135866 \cos 7 ω t + 0.00008301 \sin 7 ω t + 0.00002371 \cos 9 ω t + 0.00009085 \sin 9 ω t,

(24b)

where

ω = 1.90070426 .

(24c)

Now consider the case of β₁≠β₂. Taking each parameter as α₁ = α₂ = 1, β₁ = 1, β₂ = 999/1000, γ₁ = γ₂ = 1, and δ₁ = δ₂ = 1; in accordance with the method of the previous example, put them into (18a), (18b), (18c), and (18d), and obtain the following nonlinear algebraic equations:

6 a_{0}^{3} - 3 c_{0}^{3} + a_{0}^{2} b_{0} ω_{0} + b_{0} (- 4 + b_{0}^{2}) ω_{0} + a_{0} (4 + 6 b_{0}^{2} - 4 ω_{0}^{2}) = 0,

(25a)

6 b_{0}^{3} - a_{0} (- 4 + a_{0}^{2}) ω_{0} - a_{0} {b_{0}^{2} ω}_{0} + b_{0} (4 + 6 a_{0}^{2} - 4 ω_{0}^{2}) = 0,

(25b)

- 750 a_{0}^{3} - 750 a_{0} b_{0}^{2} + 999 c_{0} + 1500 c_{0}^{3} - 1000 c_{0} ω_{0}^{2} = 0,

(25c)

3 a_{0}^{2} b_{0} + 3 b_{0}^{3} + c_{0} (- 4 + c_{0}^{2}) ω_{0} = 0 .

(25d)

Solving the equations, obtained a₀ = 2, b₀ = 0, c₀ = 2, and ω₀ = 2.

Using the relationship between the frequency ω and the parameters ℏ is to determine the range of ℏ. Figure 5 is the curve of ω-ℏ by the third-order approximation; as can be seen from the figure, when − 0.28 < ℏ < − 0.18, the series of ω is convergent. For simplicity, take ℏ = − 0.1.

Figure 5:

The third-order approximation curve of ω versus ℏ.

Since β₁≠β₂, based on the idea of the previous example, we focus on studying the relationship between ω, a, c and the parameters β₁, β₂. First, set β₁ as an unknown parameter, β₂ as a known amount (β₂ = 999/1000); in accordance with the HAM the second-order approximate expression of ω(β₁), a(β₁), and c(β₁) can be written as

ω (β_{1}) = 2 + \frac{1}{320000000 (117997 + 25000 β_{1} + 1000 β_{1}^{2})} \times (3000000000 β_{1}^{3} + 44982001 π β_{1}^{2} - 237995998000 β_{1}^{2} + 341802016 π β_{1} - 1683662968000 β_{1} - 387215972 π + 48498065000),

(26a)

a (β_{1}) = 2 - \frac{1}{40000000 (11797 + 25000 β_{1} + 1000 β_{1}^{2})} \times (2000000000 β_{1}^{3} + 6000000000 β_{1}^{2} - 71973001 π β_{1} + 24004998000 β_{1} + 72044998 π + 112006996000),

(26b)

c (β_{1}) = 2 + \frac{1}{40000 (11797 + 25000 β_{1} + 1000 β_{1}^{2})} \times (2999 π β_{1}^{2} + 5000000 β_{1}^{2} + 65975 π β_{1} - 22000000 β_{1} - 69046 π - 126982006) .

(26c)

Next, set β₂ as an unknown parameter, β₁ as an known quantity (β₁ = 1); the second-order approximate expression of ω(β₂), a(β₂), and c(β₂) can be, respectively, written as

ω (β_{2}) = 2 - \frac{1}{320 (47 + β_{2})} \times (15 π β_{2}^{3} + 99 π β_{2}^{2} - 33 π β_{2}^{2} - 243 π β_{2} + 679 β_{2} + 129 π - 22),

(27a)

a (β_{2}) = 2 - \frac{1}{40 (47 + β_{2})} \times (π β_{2}^{3} + 21 π β_{2}^{2} - 2 β_{2}^{2} - 45 π β_{2} + 23 π + 50),

(27b)

c (β_{2}) = 2 - \frac{1}{20 (47 + β_{2})} \times (12 π β_{2}^{2} + β_{2}^{2} - 24 π β_{2} + β_{2} + 12 π + 22) .

(27c)

Draw the image of ω(β₁) and ω(β₂) in the same Cartesian coordinate system; we get curve ω-β₁ and curve ω-β₂, as Figure 6 demonstrates.

Figure 6:

The first-order approximation curves of ω-β₁ and ω-β₂.

Similar to the last case, if β₁, β₂ are within the appropriate range 1/2 < β₁, β₂ < 3/2 and the gap between β₁ and β₂ is very narrow (|β₁ − β₂| < 1/5000), then the frequency ω approximation can be stabilized at 1.95 nearby.

Figures 7 and 8 give the third-order analytical approximate phase diagram, respectively. It is not difficult to find that third-order analytical approximation agrees well with the solution obtained by the numerical integration method. Of course, to get a more accurate solution, we can appropriately increase order of count, in order to increase the number of items series solution; here the initial conditions that numerical integration method meets are x₁(0) = 2, x₁′(0) = 0, x₂(0) = 2, and x₂′(0) = 0.

Figure 7:

Comparison of the phase portrait curves of the third-order approximation with the numerical integration method.

Figure 8:

Comparison of the time history responses of the third-order approximation with the numerical integration method.

Finally, the third-order analytical approximation (x₁(t),x₂(t),ω) that the systems (1a) and (1b) satisfy the above initial conditions is

x_{1} (t) = 1.89282223 \cos ω t + 0.27163323 \sin ω t + 0.06266590 \cos 3 ω t - 0.08378610 \sin 3 ω t - 0.00331194 \cos 5 ω t - 0.00394322 \sin 5 ω t - 0.00010233 \cos 7 ω t + 0.00006081 \sin 7 ω t,

(28a)

x_{2} (t) = 1.89268859 \cos ω t + 0.27164498 \sin ω t + 0.06274216 \cos 3 ω t - 0.08378165 \sin 3 ω t - 0.00328250 \cos 5 ω t - 0.00394712 \sin 5 ω t - 0.00008900 \cos 7 ω t + 0.00005956 \sin 7 ω t,

(28b)

where

ω = 1.92252889 .

(28c)

Consider the case of γ₁≠γ₂; take α₁ = α₂ = 1, β₁ = β₂ = 1, γ₁ = 2, γ₂ = 3, and δ₁ = δ₂ = 2 and then put them into (18a), (18b), (18c), and (18d); we obtained the following nonlinear algebraic equations:

12 a_{0}^{3} - 6 c_{0}^{3} + a_{0}^{2} b_{0} ω_{0} + b_{0} (- 4 + b_{0}^{2}) ω_{0} + 4 a_{0} (1 + 3 b_{0}^{2} - 4 ω_{0}^{2}) = 0,

(29a)

4 b_{0} (1 + 3 a_{0}^{2} - 3 b_{0}^{2}) - a_{0} (- 4 + a_{0}^{2} + b_{0}^{2}) ω_{0} - 4 b_{0} ω_{0}^{2} = 0,

(29b)

9 a_{0}^{3} + 9 a_{0} b_{0}^{2} - 4 c_{0} - 15 c_{0}^{3} - 4 c_{0} ω_{0}^{2} = 0,

(29c)

9 a_{0}^{2} b_{0} + 9 b_{0}^{3} + c_{0} (- 4 + c_{0}^{2}) ω_{0} = 0 .

(29d)

Solving these equations, we obtained a₀ = 2, b₀ = 0, c₀ = 2, and $ω_{0} = \sqrt{7}$ .

Utilizing the relationship between frequency ω and parameter ℏ again is to determine the value range of ℏ. Figure 9 is the curve of ω-ℏ by the fourth-order approximation; obviously, when − 0.16 < ℏ < − 0.09, the series of ω is convergent. We may as well select ℏ = − 0.1.

Figure 9:

The fourth-order approximation curve of ω versus ℏ.

In this case, the unequal parameters γ₁, γ₂ are in the coupling, so the relationship between frequency, amplitude, and them becomes relatively simple. Set γ₁, γ₂ as unknown parameters, in accordance with the HAM; we obtained the fourth-order approximation expression of ω, a, and c. Consider the following:

ω (γ_{1}, γ_{2}) \equiv \frac{35727756872569 \sqrt{7}}{37764464640000},

(30a)

a (γ_{1}, γ_{2}) \equiv \frac{2133082117}{1105920000},

(30b)

c (γ_{1}, γ_{2}) \equiv \frac{2133082117}{1105920000} .

(30c)

Namely, as long as the order is fixed, regardless of the value of γ₁, γ₂, the approximation of frequency and amplitude is the same.

Figures 10 and 11 give the fourth-order analytical approximate phase diagram and time history curve, respectively. As is shown in the figures, fourth-order analytical approximation agrees well with the solution obtained by the numerical integration method. Here the initial conditions that numerical integration method meets are x₁(0) = 2, x₁′(0) = 0, x₂(0) = 2, and x₂′(0) = 0.

Figure 10:

Comparison of the phase portrait curves of the fourth-order approximation with the numerical integration method.

Figure 11:

Comparison of the time history responses of the fourth-order approximation with the numerical integration method.

The fourth-order analytical approximation (x₁(t),x₂(t),ω) that systems (1a) and (1b) satisfy the above initial conditions is

x_{1} (t) = 1.84397604 \cos ω t + 0.24306765 \sin ω t + 0.08670996 \cos 3 ω t - 0.06706854 \sin 3 ω t - 0.00132301 \cos 5 ω t - 0.00826489 \sin 5 ω t - 0.00057043 \cos 7 ω t - 0.00010677 \sin 7 ω t 0.00000736 \cos 9 ω t + 0.00002331 \sin 9 ω t,

(31a)

x_{2} (t) = 1.84397604 \cos ω t + 0.24306765 \sin ω t + 0.08670996 \cos 3 ω t - 0.06706854 \sin 3 ω t - 0.00132301 \cos 5 ω t - 0.00826489 \sin 5 ω t - 0.00057043 \cos 7 ω t - 0.00010677 \sin 7 ω t - 0.00000736 \cos 9 ω t + 0.00002331 \sin 9 ω t,

(31b)

ω = 2.50306103 .

(31c)

4. Four Important Theorems

From the above discussion, when α₁ = α₂ = 1, β₁ = β₂ = 1, and δ₁ = δ₂ = 2, the solution of systems (1a) and (1b) satisfies x₁(t) = x₂(t), and the expressions of x₁(t), x₂(t), and ω will not be influenced by parameters γ₁, γ₂, inspired from the results.

Theorem 1. Assuming α, β, δ, and γ_i (i = 1, 2) are real parameters, x(t) is the solution of van der Pol-Duffing equation $\ddot{x} + α (x^{2} - 1) \dot{x} + β x + δ x^{3} = 0$ ; then for any γ₁,γ₂ ∈ R, x₁(t) = x₂(t) = x(t) is a solution of system

{\ddot{x}}_{1} + α (x_{1}^{2} - 1) {\dot{x}}_{1} + β x_{1} + δ x_{1}^{3} + γ_{1} (x_{1}^{3} - x_{2}^{3}) = 0,

(32a)

{\ddot{x}}_{2} + α (x_{2}^{2} - 1) {\dot{x}}_{2} + β x_{2} + δ x_{2}^{3} + γ_{2} (x_{2}^{3} - x_{1}^{3}) = 0 .

(32b)

More generally, we have this theorem.

Theorem 2. Assuming α, β, and δ are real parameters, f_i (i = 1, 2) is real function about x₁ and x₂, and when x₁ = x₂, f₁ = f₂ = 0, x(t) is the solution of van der Pol-Duffing equation $\ddot{x} + α (x^{2} - 1) \dot{x} + β x + δ x^{3} = 0$ and x₁(t) = x₂(t) = x(t) is a solution of system

{\ddot{x}}_{1} + α (x_{1}^{2} - 1) {\dot{x}}_{1} + β x_{1} + δ x_{1}^{3} + f_{1} (x_{1}, x_{2}) = 0,

(33a)

{\ddot{x}}_{2} + α (x_{2}^{2} - 1) {\dot{x}}_{2} + β x_{2} + δ x_{2}^{3} + f_{2} (x_{1}, x_{2}) = 0 .

(33b)

It is easy to prove the above two theorems. Take x₁(t) = x₂(t) = x(t) into the systems (32a), (32b), (33a), and (33b) and the equation is established. Note that we do not have to find all solutions of the systems (32a), (32b), (33a), and (33b). However, these two theorems are of positive significance. This type of MDOF coupled van der Pol-Duffing oscillator can be simplified into a single degree of freedom of the van der Pol-Duffing equation, and then we can obtain a periodic solution.

This thought of translating MDOF into SDOF is applicable in many nonlinear systems. Here we explore another type of MDOF coupled van der Pol-Duffing oscillator by Theorems 3 and 4.

Theorem 3. Assuming α, β, δ, and γ_i (i = 1, 2) are real parameters, x(t) is the solution of van der Pol-Duffing equation $\ddot{x} + α (x^{2} - 1) \dot{x} + β x + δ x^{3} = 0$ , for any γ₁,γ₂ ∈ R, x₁(t) = x(t), and x₂(t) = − x(t) is a solution of system

{\ddot{x}}_{1} + α (x_{1}^{2} - 1) {\dot{x}}_{1} + β x_{1} + δ x_{1}^{3} + γ_{1} (x_{1}^{3} + x_{2}^{3}) = 0,

(34a)

{\ddot{x}}_{2} + α (x_{2}^{2} - 1) {\dot{x}}_{2} + β x_{2} + δ x_{2}^{3} + γ_{2} (x_{1}^{3} + x_{2}^{3}) = 0 .

(34b)

As the promotion of Theorem 3, we have Theorem 4.

Theorem 4. Assuming α, β, and δ are real parameters, f_i (i = 1, 2) is real function about x₁ and x₂, and when x₁ = x₂, f₁ = f₂ = 0, x(t) is the solution of van der Pol-Duffing equation $\ddot{x} + α (x^{2} - 1) \dot{x} + β x + δ x^{3} = 0$ , x₁(t) = x(t), and x₂(t) = − x(t) is a solution of system

{\ddot{x}}_{1} + α (x_{1}^{2} - 1) {\dot{x}}_{1} + β x_{1} + δ x_{1}^{3} + f_{1} (x_{1}, x_{2}) = 0,

(35a)

{\ddot{x}}_{2} + α (x_{2}^{2} - 1) {\dot{x}}_{2} + β x_{2} + δ x_{2}^{3} + f_{2} (x_{1}, x_{2}) = 0 .

(35b)

The proofs of Theorems 1 and 2 are exactly the same. Next we combine the HAM to study the application of Theorem 3.

In the system of (34a) and (34b), taking α = 1, β = 1, γ₁ = 5, γ₂ = 6, and δ₁ = δ₂ = 1, by Theorem 3, a periodic solution of the system (x₁(t),x₂(t)) meets x₁(t) = − x₂(t) = x(t), wherein x(t) is the solution of equation $\ddot{x} + (x^{2} - 1) \dot{x} + x + x^{3} = 0$ . Periodic solution of the equation can be approximated by using HAM. Their thoughts are similar to the previous 2-DOF coupled oscillators solving. The proof is easier, so it is not launched.

We get that the sixth-order homotopy approximate solution (x(t),ω) of the equation $\ddot{x} + (x^{2} - 1) \dot{x} + x + x^{3} = 0$ which meets the initial conditions x(0) = 2 and x′(0) = 0 is

x (t) = 1.84645601 \cos ω t + 0.33006051 \sin ω t + 0.09685555 \cos 3 ω t - 0.08775593 \sin 3 ω t - 0.00340889 \cos 5 ω t - 0.01327148 \sin 5 ω t - 0.00156607 \cos 7 ω t - 0.00028437 \sin 7 ω t - 0.00009081 \cos 9 ω t + 0.00016068 \sin 9 ω t + 0.00001575 \cos 11 ω t + 0.00001118 \sin 11 ω t + 0.00000033 \cos 13 ω t - 0.00000107 \sin 13 ω t,

(36a)

ω = 1.89931855 .

(36b)

The sixth-order homotopy approximate solution (x₁(t),x₂(t),ω) of Theorem 3 which meets the initial conditions x₁(0) = 2, x₁′(0) = 0, x₂(0) = − 2, and x₂′(0) = 0 is

x_{1} (t) = 1.84645601 \cos ω t + 0.33006051 \sin ω t + 0.09685555 \cos 3 ω t - 0.08775593 \sin 3 ω t - 0.00340889 \cos 5 ω t - 0.01327148 \sin 5 ω t - 0.00156607 \cos 7 ω t - 0.00028437 \sin 7 ω t - 0.00009081 \cos 9 ω t + 0.00016068 \sin 9 ω t + 0.00001575 \cos 11 ω t + 0.00001118 \sin 11 ω t + 0.00000033 \cos 13 ω t - 0.00000107 \sin 13 ω t,

(37a)

x_{2} (t) = - 1.84645601 \cos ω t - 0.33006051 \sin ω t - 0.09685555 \cos 3 ω t + 0.08775593 \sin 3 ω t + 0.00340889 \cos 5 ω t + 0.01327148 \sin 5 ω t + 0.00156607 \cos 7 ω t + 0.00028437 \sin 7 ω t + 0.00009081 \cos 9 ω t - 0.00016068 \sin 9 ω t - 0.00001575 \cos 11 ω t - 0.00001118 \sin 11 ω t - 0.00000033 \cos 13 ω t + 0.00000107 \sin 13 ω t,

(37b)

ω = 1.89931855 .

(37c)

Finally, we establish a comparison between the approximate solutions and the Runge-Kutta solution. Figures 12 and 13 represent a sixth-order approximation of the phase diagram and the time history graph, respectively. Clearly, the homotopy approximate solution is consistent with the Runge-Kutta solution.

Figure 12:

Comparison of the phase portrait curves of the sixth-order approximation with the numerical integration method.

Figure 13:

Comparison of the time history responses of the sixth-order approximation with the numerical integration method.

5. Conclusions

In this paper, we use the HAM to get the approximate solution of a class of coupled van der Pol oscillator. We focused on α₁≠α₂, β₁≠β₂, and γ₁≠γ₂. These three cases discussed 2-DOF coupled van der Pol system periodic solution's approximation problem, fully demonstrating the great potential of the HAM in the kinetic study. HAM is a powerful tool for solving the MDOF nonlinear system, which provides simple and ingenious ways to ensure the convergence of the series solution. It is also effective in dealing with equations without depending on the presence of a small parameter. This method not only can be obtained in any order frequency ω and expressions x₁(t) and x₂(t) but also greatly reduces the amount of the theoretical calculation. By the comparison between homotopy approximate solution and the Runge-Kutta solution, we see that the HAM can effectively deal with such parameters of coupled oscillators, just request the values of some parameters are not too big. To some degree, the proved four theorems simplify the solution of the two types of nonlinear systems. We will consider how to solve the response problem to this kind of system when the gap of parameters is relatively big in the future.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Footnotes

Acknowledgments

The authors gratefully acknowledge the support of the National Natural Science Foundations of China (NNSFC) through Grant no. 11202189 and the National Natural Science Foundation of Zhejiang through Grant no. LY12A02002.

References

Liao

S. J.

, Beyond Perturbation: Introduction to the Homotopy Analysis Method, CRC Press, Boca Raton, Fla, USA, Chapman and Hall, 2003.

Liao

S. J.

, Homotopy Analysis Method in Nonlinear Differential Equations, Springer & Higher Education Press, New York, NY, USA, 2012.

Liao

S. J.

and Chwang

A. T.

, “Application of homotopy analysis method in nonlinear oscillations,” ASME Journal of Applied Mechanics, vol. 65, no. 4, pp. 914–922, 1998.

Liao

S. J.

, “An analytic approximate approach for free oscillations of self-excited systems,” International Journal of Non-Linear Mechanics, vol. 39, no. 2, pp. 271–280, 2004.

Liao

, “On the homotopy analysis method for nonlinear problems,” Applied Mathematics and Computation, vol. 147, no. 2, pp. 499–513, 2004.

Qian

and Chen

, “Accurate approximate analytical solutions for multi-degree-of-freedom coupled van der Pol-Duffing oscillators by homotopy analysis method,” Communications in Nonlinear Science and Numerical Simulation, vol. 15, no. 10, pp. 3113–3130, 2010.

Nohara

B. T.

, and Liao

, “Series solutions of coupled Van der Pol equation by means of homotopy analysis method,” Journal of Mathematical Physics, vol. 51, no. 6, Article ID 063517, 2010.

Nayfeh

A. H.

and Mook

D. T.

, Nonlinear Oscillations, John Wiley & Sons, New York, NY, USA, 1979.

Guckenheimer

and Holmes

P. J.

, Nonlinear Oscillators, Dynamical Systems and Bifurcations of Vector Fields, Springer, New York, NY, USA, 1983.

10.

Odani

, “The limit cycle of the van der Pol equation is not algebraic,” Journal of Differential Equations, vol. 115, no. 1, pp. 146–152, 1995.

11.

Wirkus

and Rand

, “The dynamics of two coupled van der Pol oscillators with delay coupling,” Nonlinear Dynamics, vol. 30, no. 3, pp. 205–221, 2002.

12.

Hirano

and Rybicki

, “Existence of limit cycles for coupled van der Pol equations,” Journal of Differential Equations, vol. 195, no. 1, pp. 194–209, 2003.

13.

, “Dynamical analysis of two coupled parametrically excited van der Pol oscillators,” International Journal of Non-Linear Mechanics, vol. 39, no. 1, pp. 33–54, 2004.

14.

Rompala

Rand

, and Howland

, “Dynamics of three coupled van der Pol oscillators with application to circadian rhythms,” Communications in Nonlinear Science and Numerical Simulation, vol. 12, no. 5, pp. 794–803, 2007.

15.

Chen

Y. M.

and Liu

J. K.

, “A study of homotopy analysis method for limit cycle of van der Pol equation,” Communications in Nonlinear Science and Numerical Simulation, vol. 14, no. 5, pp. 1816–1821, 2009.

16.

Chen

Y. M.

and Liu

J. K.

, “Uniformly valid solution of limit cycle of the Duffing-van der Pol equation,” Mechanics Research Communications, vol. 36, no. 7, pp. 845–850, 2009.

17.

Kuznetsov

A. P.

Stankevich

N. V.

, and Turukina

L. V.

, “Coupled van der Pol-Duffing oscillators: phase dynamics and structure of synchronization tongues,” Physica D, vol. 238, no. 14, pp. 1203–1215, 2009.

18.

Kimiaeifar

Saidi

A. R.

Bagheri

G. H.

Rahimpour

, and Domairry

D. G.

, “Analytical solution for Van der Pol-Duffing oscillators,” Chaos, Solitons and Fractals, vol. 42, no. 5, pp. 2660–2666, 2009.

19.

Qian

Y. H.

Duan

C. M.

Chen

S. M.

, and Chen

S. P.

, “Asymptotic analytical solutions of the two-degree-of-freedom strongly nonlinear van der Pol oscillators with cubic couple terms using extended homotopy analysis method,” Acta Mechanica, vol. 223, no. 2, pp. 237–255, 2012.

20.

Liao

, “Notes on the homotopy analysis method: some definitions and theorems,” Communications in Nonlinear Science and Numerical Simulation, vol. 14, no. 4, pp. 983–997, 2009.