Periodic solutions of delay nonlinear system by multi-frequency homotopy analysis method

Abstract

A two-degree-of-freedom coupled Duffing system with time delay is studied by a modified homotopy analysis method. First, a detailed calculation is given for the system. Second, single periodic and period-doubling solutions are obtained by solving the constructed nonlinear algebraic equation. Third, upon comparison of the periodic solution obtained by the multi-frequency homotopy analysis method with the numerical solutions obtained by the numerical method, it is found that the approximate solution agrees well with the numerical solution for the system. Finally, we discuss the convergence control parameters and convergence control function.

Introduction

Nonlinear dynamical systems with time delay have attracted the attention of many researchers in the fields of machine tool dynamics, neural networks, biology, medicine, and population dynamics.^1–7 The time delay involved in a nonlinear dynamical system is very short in many applications, but it is still a significant factor in obtaining the correct solutions for a nonlinear system. Many researchers have pointed out that the analysis results of a nonlinear system are not correct if some small time-delay factors are ignored. For example, the time delay could affect the stability of the system, accuracy of the solution, and type of solution. Therefore, it is necessary to consider the time delay of some nonlinear dynamical systems to obtain a suitable and reasonable explanation for mechanical phenomena.

There are many ways to study time-delay systems. These methods are mainly used to investigate the periodic solution (including the single periodic solution, period-doubling solution, and quasi-periodic solution), chaos, and bifurcation behavior of a time-delay nonlinear system. Hu and Wang ⁸ described different methods of researching the dynamic behavior of time-delay systems, including the multi-scale method, method of averaging, energy analysis, and singularly perturbed pseudo-oscillator analysis. The results revealed that those methods lead to easier calculations and more accurate predictions of the local dynamics for the time-delay systems near the Hopf bifurcation. Shakeri and Dehghan ⁹ introduced the solution of differential delay equations through the homotopy perturbation method. Leung et al. ¹⁰ applied a combination of the residual harmonic balance method and polynomial homotopy continuum to analyze the steady-state bifurcation of a system with time-delay feedback control. Wang et al. ¹¹ introduced a delayed fractional financial model to obtain the effect of the delay factor on the chaos of a system, which could enhance or suppress the appearance of disorder in a system. Raghothama and Narayanan ¹² applied the incremental harmonic balance method to analyze the periodic motions of nonlinear system with time delay. Wang ¹³ presented an iterative method for calculating periodic solutions of a time-delay system caused by the Hopf bifurcation, which included degenerate cases (in which time delay disappears). Yang et al. ¹⁴ demonstrated the existence of positive periodic solutions for time-delay systems and reciprocally disturbed predator–prey systems, and gave sufficient conditions for the existence of positive periodic solutions. Su and Li ¹⁵ obtained the periodic solutions of a nonlinear parabolic system with time delay by the iterative method and the method of upper and lower solutions. Zhou and Fu ¹⁶ proved the existence of periodic solutions for a semi-linear parabolic system with discrete delay and analyzed the dynamical characteristic of the system by the method of upper and lower solutions.

The homotopy analysis method (HAM) is one of the semi-accurate methods that was presented by Liao. ¹⁷ The HAM is completely independent of small physical parameters of a system, which provides a good way to solve the family of basic functions with multiple degrees of freedom (MDOF), especially for nonlinear systems. Qian and Chen ¹⁸ applied the HAM to establish an approximate analytical solution of a nonlinear coupled oscillator with MDOF. Yuan and Li ¹⁹ calculated the first-order resonance approximation solution of the forced Duffing equation by the HAM. Shukla et al. ²⁰ obtained the periodic solutions and quasi-periodic solutions of a forced van der Pol Duffing oscillator and Duffing oscillator by the HAM. Fallahzadeh and Shakibi ²¹ acquired the series solution of the linear convection diffusion (CD) equation by the HAM. The convergence theorem is used to approximate the series solution obtained by the HAM to the exact solution of the equation. In addition, three examples were given to illustrate the efficiency and applicability of the HAM in solving the CD equation. Tajaddodianfar et al. ²² applied the HAM to obtain the first- and second-order analytical solutions for the frequency responses of the Duffing equation. Qian and Zhang ²³ introduced the optimal extended HAM to solve the damped Duffing resonator driven by a van der Pol oscillator. Results showed that the method was effective in studying the dynamic response of a MDOF nonlinear coupling system. Zou and Nagarajaiah ²⁴ proposed the multi-frequency HAM (MFHAM) to predict symmetry-breaking bifurcation and period-doubling bifurcation points accurately. Fu and Qian ²⁵ used the modified HAM to study the single periodic and period-doubling solutions of a two-DOF (TDOF) coupled Duffing system. Shen et al. ²⁶ investigated the bifurcation and chaos behaviors of a Duffing oscillator with delayed displacement and velocity feedback under harmonic excitation. Han et al. ²⁷ presented a general method for analyzing mixed-mode oscillations in parametrically and externally excited systems with two low excitation frequencies for the case of an arbitrary m:n relation between the slow frequencies of excitations. Han et al. ²⁸ also investigated the bursting dynamics of a Duffing system with multi-frequency external forcing and studied the influence of forcing amplitudes on transition of the bursting. El-Dib ²⁹ suggested a modified homotopy perturbation method by absorbing the multiple scales method. And periodical structure of the approximate solution was constructed. El-Dib ³⁰ studied the complicated dynamic problem by a coupling between the method of the multiple scales and the homotopy perturbation. Also, a second-order approximate solution was achieved.

The purpose of this paper is to investigate the periodic solution of a TDOF coupled Duffing oscillator with time-delay by the MFHAM. In addition, the dynamic response behavior of the system is analyzed.

Application of MFHAM

Consider a TDOF coupled nonlinear system with time delay

$${x^{″}}_1+c{x^{\prime }}_1+F_{1\text{τ}}(x_1,x_2)=f_1\cos \Omega t$$

(1a)

$${x^{″}}_2+c{x^{\prime }}_2+F_{2\text{τ}}(x_1,x_2)=f_2\cos \Omega t$$

(1b)where $x_i(t)(i=1,2)$ represents unknown real functions, $f_1$ and $f_2$ are the amplitudes of the parametric excitation, $F_{1\text{τ}}(x_1,x_2)$ and $F_{2\text{τ}}(x_1,x_2)$ represent coupling functions, $\Omega$ is the frequency of the parametric excitation, $c$ is a physical parameter, and $\text{τ}$ is the time delay.

Based on the MFHAM, we establish the steps of the MFHAM for the TDOF coupled nonlinear system.

An auxiliary linear-differential operator is constructed:

$$\begin{array}{c}{L}_n\left(\begin{array}{l}{x}_1\\ x_2\end{array}\right)=\left(\begin{array}{l}{x}_1^{(2n)}\\ x_2^{(2n)}\end{array}\right)+{\displaystyle \sum _i{\Omega }_i^2}\left(\begin{array}{l}{x}_1^{(2n-2)}\\ x_2^{(2n-2)}\end{array}\right)+{\displaystyle \sum _{i\ne j}{\Omega }_i^2{\Omega }_j^2}\left(\begin{array}{l}{x}_1^{(2n-4)}\\ x_2^{(2n-4)}\end{array}\right)\\ +{\displaystyle \sum _{i\ne j\ne r}{\Omega }_i^2{\Omega }_j^2{\Omega }_r^2}\left(\begin{array}{l}{x}_1^{(2n-6)}\\ x_2^{(2n-6)}\end{array}\right)+\cdots +{\displaystyle \prod _i{\Omega }_i^2}\left(\begin{array}{l}{x}_1\\ x_2\end{array}\right)\end{array}$$

(2)where ${\Omega }_i(i=1,2,\cdots )$ denotes the fundamental frequencies.

The relevant characteristic polynomial of equation (2) is

$$p(\text{λ})={\displaystyle \prod _{r=1}^n(\text{λ}+r\Omega i)}(\text{λ}-r\Omega i)$$

(3)

When the solutions of the nonlinear equation include a constant term, the auxiliary linear-differential operator (equation (2)) should be changed to:

$$\begin{array}{c}{L}_{n+1}\left(\begin{array}{l}{x}_1\\ x_2\end{array}\right)=\left(\begin{array}{l}{x}_1^{(2n+1)}\\ x_2^{(2n+1)}\end{array}\right)+{\displaystyle \sum _i{\Omega }_i^2}\left(\begin{array}{l}{x}_1^{(2n-1)}\\ x_2^{(2n-1)}\end{array}\right)+{\displaystyle \sum _{i\ne j}{\Omega }_i^2{\Omega }_j^2}\left(\begin{array}{l}{x}_1^{(2n-3)}\\ x_2^{(2n-3)}\end{array}\right)\\ +{\displaystyle \sum _{i\ne j\ne r}{\Omega }_i^2{\Omega }_j^2{\Omega }_r^2}\left(\begin{array}{l}{x}_1^{(2n-5)}\\ x_2^{(2n-5)}\end{array}\right)+\cdots +{\displaystyle \prod _i{\Omega }_i^2}\left(\begin{array}{l}{x^{\prime }}_1\\ {x^{\prime }}_2\end{array}\right)\end{array}$$

(4)

We construct the following homotopy to obtain the solution of equation (4):

$$(1-q)L_{n+1}\left(\begin{array}{l}{x}_1-g_{1,0}(t)\\ x_2-g_{2,0}(t)\end{array}\right)-q\mathit{\hslash }H(t)N\left(\begin{array}{l}{x}_1\\ x_2\end{array}\right)=0$$

(5)where $q\in [0,1]$ represents the embedding variable, and $g_{i,0}(t)(i=1,2)$ are the initial solutions of $x_1(t)$ and $x_2(t)$ , respectively. $\mathit{\hslash }$ represents the auxiliary parameter and $H(t)$ is the auxiliary function.

We assume that the power-series solutions of equation(5) are in the following form:

$$x_i=x_{i.0}+q{x}_{i,1}+q^2{x}_{i,2}+\cdots (i=1,2)$$

(6)

Substituting equation (6) into equation (5), and merging the same terms of $q$ , we obtain the following equations:

$$q^0:L_{n+1}\left(\begin{array}{l}{x}_{1,0}-g_{1,0}(t)\\ x_{2,0}-g_{2,0}(t)\end{array}\right)=0$$

(7a)

$$q^1:L_{`n+1}\left(\begin{array}{l}{x}_{1,1}-g_{1,0}(t)\\ x_{2,1}-g_{2,0}(t)\end{array}\right)=L_{n+1}\left(\begin{array}{l}{x}_{1,0}-g_{1,0}(t)\\ x_{2,0}-g_{2,0}(t)\end{array}\right)+\mathit{\hslash }H(t)N\left(\begin{array}{l}{x}_{1,0}\\ x_{2,0}\end{array}\right)$$

(7b)

$$q^2:L_{`n+1}\left(\begin{array}{l}{x}_{1,2}-g_{1,0}(t)\\ x_{2,2}-g_{2,0}(t)\end{array}\right)=L_{n+1}\left(\begin{array}{l}{x}_{1,1}-g_{1,0}(t)\\ x_{2,1}-g_{2,0}(t)\end{array}\right)+\mathit{\hslash }H(t)\frac{\partial N\left(\begin{array}{l}{x}_1\\ x_2\end{array}\right)}{\partial q}{|}_{q=0}$$

(7c)

Let the solution of equation (7a) be:

$$x_{i,0}(t)=g_{i,0}(t)+A_0+{\displaystyle \sum _{l=1}^n{A}_l\sin ({\Omega }_{l}t+{\text{φ}}_l)}$$

(8)where $A_0$ , $A_l$ , and ${\text{φ}}_l(l=1,2,\cdots ,n)$ are constants.

At this point, all of the preparation work have been finished. Now, three steps are divided to obtain the solutions of equation (5). First, substituting equation (8) into equation (7b) and eliminating the secular terms, we can obtain $L_{n+1}\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)$ . Then, substituting $L_{n+1}\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)$ into equations (2) and (4), the expressions of $x_1\left(t\right)$ and $x_2(t)$ can be obtained. Finally, we solve equations (26) and (27) including the secular terms to determine the values of $A_0$ , $A_l$ , and ${\text{φ}}_l(l=1,2,\cdots ,n)$ .

Therefore, the solutions of equation (5) are derived in the following form:

$$x_i(t)=\underset{q\to 1}{\lim }{x}_i(t,q)=x_{i,0}(t)+x_{i,1}(t)+x_{i,2}(t)+\cdots$$

(9)

Analysis of Duffing oscillator with MFHAM

Here, an example is given to illustrate the above process. A coupled Duffing system with time delay is considered:

$${x^{″}}_1(t)+c{x^{\prime }}_1(t)-x_2(t-\text{τ})+x_1^3(t)=f\cos \Omega t$$

(10a)

$${x^{″}}_2(t)+c{x^{\prime }}_2(t)-x_1(t-\text{τ})+x_2^3(t)=f\cos \Omega t$$

(10b)

Single-period solution

For the single-period solution, the characteristic polynomial is:

$$p(\text{λ})=\text{λ}{\displaystyle \prod _{r=1}^3(\text{λ}+r\Omega i)}(\text{λ}-r\Omega i)$$

(11)

The corresponding auxiliary linear-differential operator is:

$$L_4\left(\begin{array}{l}{x}_1\\ x_2\end{array}\right)=\left(\begin{array}{l}{x}_1^{(7)}\\ x_2^{(7)}\end{array}\right)+14{\Omega }^2\left(\begin{array}{l}{x}_1^{(5)}\\ x_2^{(5)}\end{array}\right)+49{\Omega }^4\left(\begin{array}{l}{x}_1^{(3)}\\ x_2^{(3)}\end{array}\right)+36{\Omega }^6\left(\begin{array}{l}{x^{\prime }}_1\\ {x^{\prime }}_2\end{array}\right)$$

(12)

The homotopy expression can be constructed in a similar way as equation (5):

$$\wp \left(\begin{array}{cc}{x}_1& q\\ x_2& q\end{array}\right)=L_4\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)-q\left(L_4\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)+\mathit{\hslash }H(t)\left(\begin{array}{c}{x^{″}}_1(t)+c{x^{\prime }}_1(t)-x_2(t-\text{τ})+x_1^3(t)-f\cos \Omega t\\ {x^{″}}_2(t)+c{x^{\prime }}_2(t)-x_1(t-\text{τ})+x_2^3(t)-f\cos \Omega t\end{array}\right)\right)$$

(13)

We assume that the solution of equation (10) can be expressed as:

$$x_i=x_{i,0}+q{x}_{i,1}+q{x}_{i,2}^2+\cdots (i=1,2)$$

(14)

Similar to equation (7), we can derive the following equations:

$$q^0:L_4\left(\begin{array}{l}{x}_{1,0}\\ x_{2,0}\end{array}\right)=0$$

(15a)

$$q^1:L_4\left(\begin{array}{l}{x}_{1,1}\\ x_{2,1}\end{array}\right)=L_4\left(\begin{array}{l}{x}_{1,0}\\ x_{2,0}\end{array}\right)+\mathit{\hslash }H(t)\left(\begin{array}{c}{x^{″}}_{1,0}(t)+c{x^{\prime }}_{1,0}(t)-x_{2,0}(t-\text{τ})+x_1^3(t)-f\cos \Omega t\\ {x^{″}}_{2,0}(t)+c{x^{\prime }}_{2,0}(t)-x_{1,0}(t-\text{τ})+x_2^3(t)-f\cos \Omega t\end{array}\right)$$

(15b)

Let the solutions of equation (15a) be expressed as

$$x_{1,0}=a_1+b_1\cos (\Omega t+{\text{φ}}_1)+b_2\cos (2\Omega t+{\text{φ}}_2)+b_3\cos (3\Omega t+{\text{φ}}_3)$$

(16a)

$$x_{2,0}=a_2+d_1\cos (\Omega t+{\text{γ}}_1)+d_2\cos (2\Omega t+{\text{γ}}_2)+d_3\cos (3\Omega t+{\text{γ}}_3)$$

(16b)where $a_1$ , $a_2$ , $b_i$ , $d_i$ , ${\text{φ}}_i$ , and ${\text{γ}}_i$ $\left(i=1,2,3\right)$ are parameters.

The expressions of $L_4(x_{i,1})(i=1,2)$ (see Appendix) can be obtained by substituting equations (16) and (15a) into equation (15b).

To eliminate the secular terms in equation (25) (see Appendix), equations (26) (see Appendix) and (27) (see Appendix) must be satisfied.

Letting $H(t)=1$ , and after elimination of the secular term, we can obtain equation (28) (see Appendix), where the parameters $a_1$ , $a_2$ , $b_i$ , $d_i$ , ${\text{φ}}_i$ , and ${\text{γ}}_i$ $\left(i=1,2,3\right)$ can be determined by equations (26) and (27).

Therefore, the first-order approximation solutions of $x_1$ and $x_2$ are:

$$x_i=\underset{q\to 1}{\lim }(x_0+q{x}_1)=x_{i,0}+x_{i,1}(i=1,2)$$

(17)

In the following investigation, the numerical simulation solutions are used to verify the solutions obtained by the MFHAM. It is found that the system has a single-periodic solution when parameters are chosen as $c=0.4,f=2.85,\Omega =1.15,\text{τ}=1.25,$ and $\mathit{\hslash }=-1$ . Figure 1 shows that the approximate solution obtained by the MFHAM agrees with the numerical solution. We then fix the parameters $f=2.88,\Omega =1.2,\text{τ}=1.25,$ and $\mathit{\hslash }=-1$ , and choose the values of $c$ as $1.5$ , $2$ and $3$ , respectively. Another three different sets of periodic solutions of the system are shown in Figures 2 to 4.

Figure 1.

Comparison of the periodic solution obtained by the MFHAM and numerical solution for $c=0.4,f=2.85,\Omega =1.15,\text{τ}=1.25,\mathit{\hslash }=-1$ . -------- numerical solution, ○○○○○○○ approximation solution. (a) Phase diagram and (b) the time history curve. MFHAM: multi-frequency homotopy analysis method.

Figure 2.

Comparison of the periodic solution obtained by the MFHAM and numerical solution for $c=1.5,f=2.88,\Omega =1.2,\text{τ}=1.25,\mathit{\hslash }=-1$ . -------- numerical solution, ○○○○○○○ approximation solution. (a) Phase diagram and (b) the time history curve.

Figure 3.

Comparison of the periodic solution obtained by the MFHAM and numerical solution for $c=2,f=2.88,\Omega =1.2,\text{τ}=1.25,\mathit{\hslash }=-1$ . -------- numerical solution, ○○○○○○○ approximation solution. (a) Phase diagram and (b) the time history curve.

Figure 4.

Comparison of the periodic solution obtained by the MFHAM and numerical solution for $c=3,f=2.88,\Omega =1.2,\text{τ}=1.25,\mathit{\hslash }=-1$ . -------- numerical solution, ○○○○○○○ approximation solution. (a) Phase diagram and (b) the time history curve.

From the above process, four different sets of parameters ( $c$ , $f$ , $\Omega$ , $\text{τ}$ and $\mathit{\hslash }$ ) are taken and four different sets of values of $a_1$ , $a_2$ , $b_i$ , $d_i$ , ${\text{φ}}_i$ , and ${\text{γ}}_i$ $(i=1,2,3)$ are obtained as shown in Tables 1 to 4, respectively.

Table 1.

The value of each parameter by eliminating the secular term, in the case of $c=0.4,f=2.85,\Omega =1.15,\tau =1.25,\mathit{\hslash }=-1$ .

$a_1$	$a_2$	$b_1$	$b_2$	$b_3$	$d_1$	$d_2$
$\begin{array}{l}-1.2764\\ \times {10}^{-18}\end{array}$	$\begin{array}{l}-1.27641\\ \times {10}^{-18}\end{array}$	$1.67089$	$\begin{array}{l}5.2054\\ \times {10}^{-20}\end{array}$	$\begin{array}{l}-0.1593\\ 06\end{array}$	$1.67089$	$\begin{array}{l}5.20633\\ \times {10}^{-20}\end{array}$
$d_3$	${\text{φ}}_1$	${\text{φ}}_2$	${\text{φ}}_3$	${\text{γ}}_1$	${\text{γ}}_2$	${\text{γ}}_3$
$\begin{array}{l}-0.1593\\ 06\end{array}$	$-1.03156$	$30.7565$	$37.8087$	$-1.03156$	$30.7565$	$37.8087$

Table 2.

The value of each parameter by eliminating the secular term, in the case of $c=1.5,f=2.88,\Omega =1.2,\tau =1.25,\mathit{\hslash }=-1$ .

$a_1$	$a_2$	$b_1$	$b_2$	$b_3$	$d_1$	$d_2$
$\begin{array}{l}-3.03024\\ \times {10}^{-25}\end{array}$	$\begin{array}{l}-3.06962\\ \times {10}^{-25}\end{array}$	$\begin{array}{l}-0.9799\\ 51\end{array}$	$\begin{array}{l}-4.32674\\ \times {10}^{-26}\end{array}$	$\begin{array}{l}-0.0193\\ 752\end{array}$	$\begin{array}{l}-0.9799\\ 51\end{array}$	$\begin{array}{l}-4.22703\\ \times {10}^{-26}\end{array}$
$d_3$	${\text{φ}}_1$	${\text{φ}}_2$	${\text{φ}}_3$	${\text{γ}}_1$	${\text{γ}}_2$	${\text{γ}}_3$
$\begin{array}{l}-0.0193\\ 752\end{array}$	$-4.98273$	$-24.6998$	$-24.0002$	$-4.98273$	$-24.6998$	$-24.0002$

Table 3.

The value of each parameter by eliminating the secular term, in the case of $c=2,f=2.88,\Omega =1.2,\tau =1.25,\mathit{\hslash }=-1$ .

$a_1$	$a_2$	$b_1$	$b_2$	$b_3$	$d_1$	$d_2$
$\begin{array}{l}7.14695\\ \times {10}^{-25}\end{array}$	$\begin{array}{l}7.14695\\ \times {10}^{-25}\end{array}$	$0.802693$	$\begin{array}{l}-1.11539\\ \times {10}^{-25}\end{array}$	$\begin{array}{l}-0.0097\\ 0348\end{array}$	$0.802693$	$\begin{array}{l}-1.11551\\ \times {10}^{-25}\end{array}$
$d_3$	${\text{φ}}_1$	${\text{φ}}_2$	${\text{φ}}_3$	${\text{γ}}_1$	${\text{γ}}_2$	${\text{γ}}_3$
$\begin{array}{l}-0.0097\\ 0348\end{array}$	$-1.86283$	$6.61507$	$-8.2442$	$-1.86283$	$6.61507$	$-8.2442$

Table 4.

The value of each parameter by eliminating the secular term, in the case of $c=3,f=2.88,\Omega =1.2,\tau =1.25,\mathit{\hslash }=-1$ .

$a_1$	$a_2$	$b_1$	$b_2$	$b_3$	$d_1$	$d_2$
$\begin{array}{l}-0.6596\\ 01\end{array}$	$\begin{array}{l}-0.6596\\ 01\end{array}$	$\begin{array}{l}0.6079\\ 34\end{array}$	$\begin{array}{l}-0.0465\\ 598\end{array}$	$\begin{array}{l}0.0062\\ 6846\end{array}$	$\begin{array}{l}0.6079\\ 34\end{array}$	$\begin{array}{l}-0.0465\\ 598\end{array}$
$d_3$	${\text{φ}}_1$	${\text{φ}}_2$	${\text{φ}}_3$	${\text{γ}}_1$	${\text{γ}}_2$	${\text{γ}}_3$
$\begin{array}{l}0.0062\\ 6846\end{array}$	$-1.54759$	$-20.7311$	$-22.142$	$-1.54759$	$-20.7311$	$-22.142$

Period-doubling solution

For the period-doubling solution of equation (10), the characteristic polynomial is:

$$p\left(\text{λ}\right)=\text{λ}{\displaystyle \prod _{r=1}^6\left(\text{λ}+\frac{r\Omega }{2}i\right)\left(\text{λ}-\frac{r\Omega }{2}i\right)}$$

(18)and the corresponding auxiliary linear-differential operator is of the following form:

$$\begin{array}{c}{L}_7\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)=\left(\begin{array}{c}{x}_1^{\left(13\right)}\\ x_2^{\left(13\right)}\end{array}\right)+\frac{91{\Omega }^2}{4}\left(\begin{array}{c}{x}_1^{\left(11\right)}\\ x_2^{\left(11\right)}\end{array}\right)+\frac{3003{\Omega }^4}{16}\left(\begin{array}{c}{x}_1^{\left(9\right)}\\ x_2^{\left(9\right)}\end{array}\right)+\frac{44473{\Omega }^6}{64}\left(\begin{array}{c}{x}_1^{\left(7\right)}\\ x_2^{\left(7\right)}\end{array}\right)\\ +\frac{37037{\Omega }^8}{32}\left(\begin{array}{c}{x}_1^{\left(5\right)}\\ x_2^{\left(5\right)}\end{array}\right)+\frac{48321{\Omega }^{10}}{64}\left(\begin{array}{c}{x}_1^{\left(3\right)}\\ x_2^{\left(3\right)}\end{array}\right)+\frac{2025{\Omega }^{10}}{16}\left(\begin{array}{c}{x^{\prime }}_1\\ {x^{\prime }}_2\end{array}\right)\end{array}$$

(19)

The homotopy expression, which is similar to that obtained in Section ‘Single-period solution’, is:

$$\wp \left(\begin{array}{cc}{x}_1& q\\ x_2& q\end{array}\right)=L_7\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)-q\left(L_7\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)+\mathit{\hslash }H\left(t\right)\left(\begin{array}{c}{x^{″}}_1(t)+c{x^{\prime }}_1(t)-x_2(t-\text{τ})+x_1^3(t)-f\cos \Omega t\\ {x^{″}}_2(t)+c{x^{\prime }}_2(t)-x_1(t-\text{τ})+x_2^3(t)-f\cos \Omega t\end{array}\right)\right)$$

(20)

We assume that the solution of equation (10) can be expressed as

$$x_i=x_{i,0}+q{x}_{i,1}+q^2{x}_{i,2}+\cdots \left(i=1,2\right)$$

(21)

Substituting equation (21) into equation (20), and merging the same term of $q$ , we can obtain:

$$q^0:L_7\left(\begin{array}{c}{x}_{1,0}\\ x_{2,0}\end{array}\right)=0$$

(22a)

$$q^1:L_7\left(\begin{array}{l}{x}_{1,1}\\ x_{2,1}\end{array}\right)=L_7\left(\begin{array}{l}{x}_{1,0}\\ x_{2,0}\end{array}\right)+\mathit{\hslash }H(t)\left(\begin{array}{c}{x^{″}}_{1,0}(t)+c{x^{\prime }}_{1,0}(t)-x_{2,0}(t-\text{τ})+x_1^3(t)-f\cos \Omega t\\ {x^{″}}_{2,0}(t)+c{x^{\prime }}_{2,0}(t)-x_{1,0}(t-\text{τ})+x_2^3(t)-f\cos \Omega t\end{array}\right)$$

(22b)

We suppose that the solutions of equation (22a) are:

$$\begin{array}{c}{x}_{1,0}=a_1+b_1\cos (\Omega t+{\text{φ}}_1)+b_2\cos (2\Omega t+{\text{φ}}_2)+b_3\cos (3\Omega t+{\text{φ}}_3)\\ +e_1\cos \left(\frac{\Omega t}{2}+{\text{φ}}_0\right)+e_2\cos \left(\frac{3\Omega t}{2}+{\text{φ}}_4\right)+e_3\cos \left(\frac{5\Omega t}{2}+{\text{φ}}_5\right)\end{array}$$

(23a)

$$\begin{array}{c}{x}_{2,0}=a_2+d_1\cos (\Omega t+{\text{φ}}_1)+d_2\cos (2\Omega t+{\text{φ}}_2)+d_3\cos (3\Omega t+{\text{φ}}_3)\\ +f_1\cos \left(\frac{\Omega t}{2}+{\text{γ}}_0\right)+f_2\cos \left(\frac{3\Omega t}{2}+{\text{γ}}_4\right)+f_3\cos \left(\frac{5\Omega t}{2}+{\text{γ}}_5\right)\end{array}$$

(23b)where $a_1,a_2,b_i,d_i,e_i,f_i(i=1,2,3),$ ${\text{φ}}_j,$ and ${\text{γ}}_j(j=0,1,2,3,4,5)$ are parameters.

By using the same step as presented in Section ‘Single-period solution’, the period-doubling solution is obtained.

Choosing $c=0.2,f=4.8,\Omega =1.2,\text{τ}=1.25,$ and $\mathit{\hslash }=-1$ , the system has the period-doubling solution. The good agreement between the approximate and numerical solutions is as shown in Figure 5. The corresponding values of parameters $a_1,a_2,b_i,d_i,e_i,f_i(i=1,2,3)\text{,}$ ${\text{φ}}_j,$ and ${\text{γ}}_j(j=0,1,2,3,4,5)$ are given in Table 5.

Figure 5.

Comparison of the period-doubling solution obtained by the MFHAM and numerical solution for $c=0.2,f=4.8,\Omega =1.2,\text{τ}=1.25,\mathit{\hslash }=-1$ . -------- numerical solution, ○○○○○○○ approximation solution. (a) phase diagram and (b) the time history curve.

Table 5.

The value of each parameter by eliminating the secular term, in the case of $c=0.2,f=4.8,\Omega =1.2,\text{τ}=1.25,\mathit{\hslash }=-1$ .

$a_1$	$a_2$	$b_1$	$b_2$	$b_3$	$d_1$	$d_2$
$0.415039$	$0.415039$	$1.81145$	$1.02152$	$\begin{array}{l}-0.0461\\ 115\end{array}$	$1.81145$	$1.02152$
$d_3$	$e_1$	$e_2$	$e_3$	$f_1$	$f_2$	$f_3$
$\begin{array}{l}-0.0461\\ 115\end{array}$	$\begin{array}{l}-1.18474\\ \times {10}^{-32}\end{array}$	$\begin{array}{l}7.70372\\ \times {10}^{-33}\end{array}$	$\begin{array}{l}-3.8441\\ 6\times {10}^{-31}\end{array}$	$\begin{array}{l}1.76522\\ \times {10}^{-31}\end{array}$	$\begin{array}{l}-1.0323\\ \times {10}^{-31}\end{array}$	$\begin{array}{l}5.93186\\ \times {10}^{-32}\end{array}$
${\text{φ}}_0$	${\text{φ}}_1$	${\text{φ}}_2$	${\text{φ}}_3$	${\text{φ}}_4$	${\text{φ}}_5$	${\text{γ}}_0$
$2.59959$	$\begin{array}{l}-0.7421\\ 29\end{array}$	$1.17113$	$0.504892$	$3.19887$	$1.09963$	$2.59959$
${\text{γ}}_1$	${\text{γ}}_2$	${\text{γ}}_3$	${\text{γ}}_4$	${\text{γ}}_5$
$\begin{array}{l}-0.7421\\ 29\end{array}$	$1.17113$	$0.504892$	$3.19887$	$1.09963$

Comments on homotopy analysis method

Here, we present a nonlinear system and discuss the convergence control of the solution by taking different convergence parameters and convergence functions.

We consider the following system:

$${x^{″}}_1(t)+c{x^{\prime }}_1(t)-x_2(t-\text{τ})+x_1^3(t)=f\cos \Omega t$$

(24a)

$${x^{″}}_2(t)+c{x^{\prime }}_2(t)-x_1(t-\text{τ})+x_2^3(t)=f\cos \Omega t$$

(24b)

First, we fix the parameters $c=0.4,f=2.85,\Omega =1.15$ , and $\tau =1.25$ , and further take $\mathit{\hslash }=0.0001,h=1,h=1000,h=10\text{,}000$ , and $\mathit{\hslash }=100,000$ , respectively, Figures 6 to 10 compare the approximate and exact solutions of equation (24). We find that the three images obtained by $\mathit{\hslash }=0.0001$ (Figure 6), $\mathit{\hslash }=1$ (Figure 7) and $\mathit{\hslash }=1000$ (Figure 8) are basically the same, and the approximate solutions of $x_1(t)$ and $x_2(t)$ are in good agreement with the exact solution. However, for $\mathit{\hslash }=10,000$ (Figure 9), the approximate solutions of $x_1(t)$ and $x_2(t)$ do not agree well with the exact solution. For $\mathit{\hslash }=100,000$ (Figure 10), the approximate solutions of $x_1(t)$ and $x_2(t)$ do not agree with the exact solution.

Figure 6.

For $\mathit{\hslash }=0.0001$ , $c=0.4,f=2.85,\Omega =1.15,\text{τ}=1.25$ , comparison of the approximate solution and exact solution for equation (24). -------- numerical solution, ○○○○○○○ approximation solution.

Figure 7.

For $\mathit{\hslash }=1$ , $c=0.4,f=2.85,\Omega =1.15,\text{τ}=1.25$ , comparison of the approximate solution and exact solution for equation (24). -------- numerical solution, ○○○○○○○ approximation solution.

Figure 8.

For $\mathit{\hslash }=1000$ , $c=0.4,f=2.85,\Omega =1.15,\text{τ}=1.25$ , comparison of the approximate solution and exact solution for equation (24). -------- numerical solution, ○○○○○○○ approximation solution.

Figure 9.

For $\mathit{\hslash }=10000$ , $c=0.4,f=2.85,\Omega =1.15,\text{τ}=1.25$ , comparison of the approximate solution and exact solution for equation (24). -------- numerical solution, ○○○○○○○ approximation solution.

Figure 10.

For $\mathit{\hslash }=100000$ , $c=0.4,f=2.85,\Omega =1.15,\text{τ}=1.25$ , comparison of the approximate solution and exact solution for equation (24). -------- numerical solution, ○○○○○○○ approximation solution.

Therefore, in summary, when $\mathit{\hslash }$ is greater than or equal to 10,000, the deviation between the approximate and exact solutions also increases gradually.

Conclusion

In this paper, we applied the MFHAM to study the nonlinear dynamic behaviors of the TDOF coupled Duffing system with time delay. Here, not only have the steps of solving the system been constructed based on the MFHAM, but it also has been found that the single-period and period-doubling solutions obtained by the MFHAM are in good agreement with the numerical solutions. Therefore, the MFHAM is an effective method used to solve the periodic solutions of the TDOF coupled Duffing oscillator with time delay. In Section ‘Comments on homotopy analysis method’, we found that when $\mathit{\hslash }$ is greater than a certain value, the approximate solution of equation (24) does not agree with the exact solution.