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

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

x ¨ 1 + α 1 ( x 1 2 - 1 ) x ˙ 1 + β 1 x 1 + δ 1 x 1 3 + γ 1 ( x 1 3 - x 2 3 ) = 0 ,

(1a)

x ¨ 2 + α 2 ( x 2 2 - 1 ) 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 ( τ ) = ∑ k = 0 + ∞ ( a 1 , k cos k τ + b 1 , k sin k τ ) ,

(4a)

x 2 ( τ ) = ∑ k = 0 + ∞ ( 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 ( Φ 1 ( τ , q ) Φ 2 ( τ , q ) ) = ω 0 2 ( ∂ 2 Φ 1 ( τ , q ) ∂ τ 2 + Φ 1 ( τ , q ) ∂ 2 Φ 2 ( τ , q ) ∂ τ 2 + Φ 2 ( τ , q ) ) .

(6)

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

N ( Φ 1 ( τ , q ) Φ 2 ( τ , q ) ) = ( Ω 2 ( q ) ∂ 2 Φ 1 ( τ , q ) ∂ τ 2 + Ω ( q ) α 1 [ Φ 1 2 ( τ , q ) - 1 ] ∂ Φ 1 ( τ , q ) ∂ τ + β 1 Φ 1 ( τ , q ) + δ 1 Φ 1 3 ( τ , q ) + γ 1 [ Φ 1 3 ( τ , q ) - Φ 2 3 ( τ , q ) ] Ω 2 ( q ) ∂ 2 Φ 2 ( τ , q ) ∂ τ 2 + Ω ( q ) α 2 [ Φ 2 2 ( τ , q ) - 1 ] ∂ Φ 2 ( τ , q ) ∂ τ + β 2 Φ 2 ( τ , q ) + δ 2 Φ 2 3 ( τ , q ) + γ 2 [ Φ 2 3 ( τ , q ) - Φ 1 3 ( τ , q ) ] ) .

(7)

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

( 1 - q ) L ( Φ 1 ( τ , q ) - x 1 , 0 ( τ ) Φ 2 ( τ , q ) - x 2 , 0 ( τ ) ) = q ℏ N ( Φ 1 ( τ , q ) Φ 2 ( τ , q ) ) ,

(8)

with the initial conditions

Φ 1 ( 0 , q ) = a ( q ) , ∂ Φ 1 ( τ , q ) ∂ τ | τ = 0 = b ( q ) ,

(9a)

Φ 2 ( 0 , q ) = c ( q ) , ∂ Φ 2 ( τ , q ) ∂ τ | τ = 0 = 0 .

(9b)

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

Φ 1 ( τ , 0 ) = x 1 , 0 ( τ ) , Φ 2 ( τ , 0 ) = x 2 , 0 ( τ ) , Ω ( 0 ) = ω 0 .

(10)

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

Φ 1 ( τ , 1 ) = x 1 ( τ ) , Φ 2 ( τ , 1 ) = x 2 ( τ ) , Ω ( 1 ) = ω .

(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 ( τ ) + ∑ m = 1 + ∞ x 1 , m ( τ ) q m ,

(12a)

Φ 2 ( τ , q ) = x 2 , 0 ( τ ) + ∑ m = 1 + ∞ x 2 , m ( τ ) q m ,

(12b)

Ω ( q ) = ω 0 + ∑ m = 1 + ∞ ω m q m ,

(12c)

x 1 , m ( τ ) = 1 m ! ∂ m Φ 1 ( τ , q ) ∂ q m | q = 0 ,

(13a)

x 2 , m ( τ ) = 1 m ! ∂ m Φ 2 ( τ , q ) ∂ q m | q = 0 ,

(13b)

ω m = 1 m ! ∂ m Ω ( q ) ∂ 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 ( τ ) + ∑ m = 1 + ∞ x 1 , m ( τ ) ,

(14a)

x 2 ( τ ) = x 2 , 0 ( τ ) + ∑ m = 1 + ∞ x 2 , m ( τ ) ,

(14b)

ω = ω 0 + ∑ m = 1 + ∞ ω m .

(14c)

For convenience, we define the vector as follows:

x 1 , n = { x 1 , 0 ( τ ) , x 1 , 1 ( τ ) , … , x 1 , n ( τ ) } ,

(15a)

x 2 , n = { x 2 , 0 ( τ ) , x 2 , 1 ( τ ) , … , x 2 , n ( τ ) } ,

(15b)

ω n = { ω 0 , ω 1 , … , ω 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 ( x 1 , m ( τ ) - χ m x 1 , m - 1 ( τ ) x 2 , m ( τ ) - χ m x 2 , m - 1 ( τ ) ) = ℏ ( R 1 , m ( x 1 , m - 1 , ω m - 1 ) R 2 , m ( x 2 , m - 1 , ω m - 1 ) ) ,

(16)

with the initial conditions

x 1 , m ( 0 ) = a m , x 1 , m ′ ( 0 ) = b m , x 2 , m ( 0 ) = c m , x 2 , m ′ ( 0 ) = 0 , ( m ≥ 1 ) , ( R 1 , m ( x 1 , m - 1 , ω m - 1 ) R 2 , m ( x 2 , m - 1 , ω m - 1 ) ) = 1 ( m - 1 ) ! ∂ m - 1 N ( Φ 1 ( τ , q ) Φ 2 ( τ , q ) ) ∂ q m - 1 | q = 0 .

(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

1 π ∫ 0 2 π R 1 , m ( x 1 , m - 1 , ω m - 1 ) cos τ d τ = 0 ,

(18a)

1 π ∫ 0 2 π R 2 , m ( x 2 , m - 1 , ω m - 1 ) cos τ d τ = 0 ,

(18b)

1 π ∫ 0 2 π R 1 , m ( x 1 , m - 1 , ω m - 1 ) sin τ d τ = 0 ,

(18c)

1 π ∫ 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 ) = ∑ 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 ) = ∑ n = 1 φ 2 ( m ) [ a m , n 2 cos n τ + b m , n 2 sin n τ ] ,

(19b)

a m , n 1 = 1 π ∫ 0 2 π R 1 , m ( x 1 , m - 1 , ω m - 1 ) cos ( n τ ) d τ ,

(20a)

a m , n 2 = 1 π ∫ 0 2 π R 2 , m ( x 2 , m - 1 , ω m - 1 ) cos ( n τ ) d τ ,

(20b)

b m , n 1 = 1 π ∫ 0 2 π R 1 , m ( x 1 , m - 1 , ω m - 1 ) sin ( n τ ) d τ ,

(20c)

b m , n 2 = 1 π ∫ 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.

OPEN IN VIEWER

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 ) = 44779 256000000 α 1 4 - 1731873 6400000000 α 1 3 - 23516728323 1280000000000 α 1 2 + 986875659 102400000000 α 1 + 78481134113267 40960000000000 ,

(22a)

a ( α 1 ) = - 121 4000000 α 1 4 + 714791 800000000 α 1 3 - 2575699901 240000000000 α 1 2 + 81116409 2560000000 α 1 + 2458329452757 1280000000000 ,

(22b)

c ( α 1 ) = 13 115200 α 1 4 - 363 1280000 α 1 3 - 3893383 384000000 α 1 2 + 138557859 12800000000 α 1 + 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 ) = 1937 11059200 α 2 4 - 757 2764800 α 2 3 - 101519 5529600 α 2 2 + 2989 307200 α 2 + 784659 409600 ,

(23a)

a ( α 2 ) = 17 207360 α 2 4 + 919 1036800 α 2 3 - 23207 1036800 α 2 2 + 3691 115200 α 2 + 1159 600 ,

(23b)

c ( α 2 ) = - 31 115200 α 2 3 + 19 12800 α 2 2 + 419 38400 α 2 + 44471 23040 .

(23c)

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

OPEN IN VIEWER

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.

OPEN IN VIEWER

Figure 3:

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

OPEN IN VIEWER

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.

OPEN IN VIEWER

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 + 1 320000000 ( 117997 + 25000 β 1 + 1000 β 1 2 ) × ( 3000000000 β 1 3 + 44982001 π β 1 2 - 237995998000 β 1 2 + 341802016 π β 1 - 1683662968000 β 1 - 387215972 π + 48498065000 ) ,

(26a)

a ( β 1 ) = 2 - 1 40000000 ( 11797 + 25000 β 1 + 1000 β 1 2 ) × ( 2000000000 β 1 3 + 6000000000 β 1 2 - 71973001 π β 1 + 24004998000 β 1 + 72044998 π + 112006996000 ) ,

(26b)

c ( β 1 ) = 2 + 1 40000 ( 11797 + 25000 β 1 + 1000 β 1 2 ) × ( 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 - 1 320 ( 47 + β 2 ) × ( 15 π β 2 3 + 99 π β 2 2 - 33 π β 2 2 - 243 π β 2 + 679 β 2 + 129 π - 22 ) ,

(27a)

a ( β 2 ) = 2 - 1 40 ( 47 + β 2 ) × ( π β 2 3 + 21 π β 2 2 - 2 β 2 2 - 45 π β 2 + 23 π + 50 ) ,

(27b)

c ( β 2 ) = 2 - 1 20 ( 47 + β 2 ) × ( 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.

OPEN IN VIEWER

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.

OPEN IN VIEWER

Figure 7:

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

OPEN IN VIEWER

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 = 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.

OPEN IN VIEWER

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 ) ≡ 35727756872569 7 37764464640000 ,

(30a)

a ( γ 1 , γ 2 ) ≡ 2133082117 1105920000 ,

(30b)

c ( γ 1 , γ 2 ) ≡ 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.

OPEN IN VIEWER

Figure 10:

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

OPEN IN VIEWER

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 x ¨ + α ( x 2 - 1 ) x ˙ + β x + δ x 3 = 0 ; then for any γ₁,γ₂ ∈ R, x₁(t) = x₂(t) = x(t) is a solution of system

x ¨ 1 + α ( x 1 2 - 1 ) x ˙ 1 + β x 1 + δ x 1 3 + γ 1 ( x 1 3 - x 2 3 ) = 0 ,

(32a)

x ¨ 2 + α ( x 2 2 - 1 ) 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 x ¨ + α ( x 2 - 1 ) x ˙ + β x + δ x 3 = 0 and x₁(t) = x₂(t) = x(t) is a solution of system

x ¨ 1 + α ( x 1 2 - 1 ) x ˙ 1 + β x 1 + δ x 1 3 + f 1 ( x 1 , x 2 ) = 0 ,

(33a)

x ¨ 2 + α ( x 2 2 - 1 ) 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 x ¨ + α ( x 2 - 1 ) x ˙ + β x + δ x 3 = 0 , for any γ₁,γ₂ ∈ R, x₁(t) = x(t), and x₂(t) = − x(t) is a solution of system

x ¨ 1 + α ( x 1 2 - 1 ) x ˙ 1 + β x 1 + δ x 1 3 + γ 1 ( x 1 3 + x 2 3 ) = 0 ,

(34a)

x ¨ 2 + α ( x 2 2 - 1 ) 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 x ¨ + α ( x 2 - 1 ) x ˙ + β x + δ x 3 = 0 , x₁(t) = x(t), and x₂(t) = − x(t) is a solution of system

x ¨ 1 + α ( x 1 2 - 1 ) x ˙ 1 + β x 1 + δ x 1 3 + f 1 ( x 1 , x 2 ) = 0 ,

(35a)

x ¨ 2 + α ( x 2 2 - 1 ) 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 x ¨ + ( x 2 - 1 ) 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 x ¨ + ( x 2 - 1 ) 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.

OPEN IN VIEWER

Figure 12:

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

OPEN IN VIEWER

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.