Triple-Zero Bifurcation of Van Der Pol Oscillator with Delay Feedback

Suqi Ma

Department of Mathematics, China Agricultural University, Beijing, China.

DOI: 10.4236/ijmnta.2021.103008                

Abstract

A van der Pol equation underlying state feedback control is investigated and the triple-zero bifurcation arises at the bifurcation point which is of codimension three singularity. By applying Schmidt-Lyapunov reduction method combined with center manifold analytical technique, the near approximating formal norm is derived at the triple-zero point. Hence after, as varying parameters continuously, the numerical simulation produces homoclinic bifurcation solutions appearing in system. In addition, the numerical simulation also exhibits the produced double-period limit cycle with chosen bifurcation parameters and the routes to chaos via period-doubling bifurcation are also verified.

Keywords

Triple-Zero Bifurcation, Homoclinic Solutions, DDEs

1. Introduction

Within the discipline of dynamical interests, the dynamics of van der Pol oscillator focused most attention since disparate sorts of complex dynamical phenomena appear in every field of scientific investigation [1] [2] . For example, it is often exemplified as a prototype equation to be applied in the investigation work of limit cycle oscillation with self-excitation [3] . Especially, people have made big effort to study the complex bifurcation behavior of van der Pol oscillator with state feedback control with the inherent discrete time delay, which may further bring complex dynamics and high-codimensional bifurcation behavior. For example, Hopf bifurcation arose with the loss of stability of the trivial solution in van der Pol oscillator and it has been reported that delay feedback may change the amplitudes of the associated limit cycle [4] [5] . The high-codimensional bifurcation may bring forth torus, chaos and strange attractors in system [6] [7] .

Triple-zero point represents that the bifurcation of the related equilibrium solution has algebra multiple-3 and geometrical multiple-1 singularity. Generally system dynamical behavior with triple-zero bifurcation may be induced due to Hopf bifurcation and double zero bifurcation merging in the triple-zero point [6] [8] . Triple-zero bifurcation in van der Pol oscillator with delay feedback has been reported in paper [9] [10] ; both Fold Hopf bifurcation and Takens-Bogdanov bifurcation are discovered near the triple-zero three singularity bifurcation point. It is pointed out, both limit cycles with double period and chaos are produced in the neighbor regime of the triple zero point [11] [12] [13] . In addition, the associated homoclinic bifurcation is always highlighted in the studying of dynamical behavior brought forth by triple-zero bifurcation and the related bifurcation phenomena. However, the similar bifurcation behavior of van der Pol oscillator with commensurable two-time delay feedbacks control is still few examples. The investigation work further explores high codimension bifurcation behavior with state feedback control acted by multiple time delay. The main feature manifest in triple-zero bifurcation is still unclear and needs more effort.

In this paper the dynamics of a general van der Pol oscillator with two commensurable time delay in state feedback control is discussed, and the corresponding system is listed as

$\ddot{x}+\epsilon \left(x^2-1\right)\stackrel{\dot{}}{x}+x=\epsilon kx\left(t-\tau \right)+\epsilon kx\left(t-2\tau \right)-\beta {\stackrel{\dot{}}{x}}^2+\gamma x^4\mathrm{,}x\in R$ (1.1)

With the background knowledge of delay differential equations, we know system (1.1) is the infinite dimensional system. It is defined on the continuous functional space $C=C^{\prime }\left(\left[-2\tau \mathrm{,0}\right]\mathrm{,}{R}^2\right)$ which equipped with the norm $\Vert \varphi \Vert ={\sup }_{\theta \in \left[-2\tau \mathrm{,0}\right]}\left|\varphi \left(\theta \right)\right|$ . System (1.1) can experience the cascades of bifurcation behaviors to chaos phenomena. With free parameters $\epsilon \mathrm{,}k\mathrm{,}\tau$ , the routes via period doubling bifurcation of period solutions to chaos are explored.

With the aids of parameter perturbation near triple-zero bifurcation point, the normal form of system (1.1) is analyzed by using dimension reduction method combined with center manifold analytical theory in DDEs [14] [15] . With the triple zero bifurcation phenomena, nearby Hopf bifurcation and double zero bifurcation are detected and analyzed in parameter space with the norm form theory [16] , which further explains the complex dynamical phenomena classified in parameter space. The derived universal unfolding of normal form is useful in understanding the local bifurcation behavior as Takens-Bogdonav singularity and homoclinic solution bifurcated from the triple-zero point. Accordingly, the bifurcated homoclinic solution of DDEs (1.1) is computed and verified by DDE-Biftool simulation software [17] [18] . As varying time delay, the double-limit cycle is produced near bifurcation point of high codimension singularity.

The whole paper is organized and listed as the following: In Section 2, the triple zero bifurcation of system (1.1) is analyzed with delay state feedback control and the triple zero singularity is discussed in Section 3. In Section 4, the normal form near the triple zero point is derived by using center manifold theory and parameter perturbation method. Henceforth, the homoclinic solutions which emanating from triple-zero point are computed by varying free parameters continuously and the routes from double period solutions to chaos are simulated by choosing Poincare section underlying time delay control.

2. Triple Zero Bifurcation of Continuous System

To discuss the triple-zero bifurcation of system (1.1), roots with zero real parts of the related characteristic equation are calculated by the linearized equation of DDEs (1.1). Equation (1.1) is rewritten as

$\begin{array}{l}\stackrel{\dot{}}{x}=y,\\ \stackrel{\dot{}}{y}=-\epsilon \left(x^2-1\right)y-x+\epsilon kx\left(t-\tau \right)+\epsilon kx\left(t-2\tau \right)-\beta y^2+\gamma x^4,\end{array}$ (2.1)

The linearized equation of Equations (2.1) is

$\begin{array}{l}\stackrel{\dot{}}{x}=y,\\ \stackrel{\dot{}}{y}=\epsilon y-x+\epsilon kx\left(t-\tau \right)+\epsilon kx\left(t-2\tau \right),\end{array}$ (2.2)

For any $\varphi \in C$ , by Rieze representation theorem, there exists $2\times 2$ bounded variation function $\eta \left(\theta \right)$ to express Equation (2.2) as linear operator

$L\varphi ={\displaystyle {\int }_{-\tau }^0}\text{d}\eta \left(\theta \right)\varphi \left(\theta \right)$ (2.3)

with

$\text{d}\eta =\left(\begin{array}{cc}0& \delta \left(\theta \right)\\ -\delta \left(\theta \right)+\epsilon k\delta \left(\theta +\tau \right)+\epsilon k\delta \left(\theta +2\tau \right)& \epsilon \delta \left(\theta \right)\end{array}\right)$

The related characteristic equation is written as

$\Delta \left(\lambda \right)=\lambda -{\displaystyle {\int }_{-2\tau }^0}\text{d}\eta \left(\theta \right)\exp \left(\lambda \theta \right)$ (2.4)

that is,

$\Delta \left(\lambda \right)={\lambda }^2-\epsilon \lambda +1-\epsilon k{\text{e}}^{-\lambda \tau }-\epsilon k{\text{e}}^{-2\lambda \tau }\equiv 0$ (2.5)

Obviously, the fold bifurcation occurs at the line $2k\epsilon =1$ and TB bifurcation occurs at point $\left(\epsilon \mathrm{,}k\mathrm{,}\tau \right)=\left(\epsilon \mathrm{,}\frac{1}{2\epsilon }\mathrm{,}\frac{2\epsilon }{3}\right)$ since a double zero eigenvalues with algebraic multiplicity of 2 appearing in the characteristic Equation (2.5). Set $\lambda =\pm i\omega \left(\omega >0\right)$ , then we have

$\begin{array}{l}-{\omega }^2+1=\epsilon k\cos \left(\omega \tau \right)+\epsilon k\cos \left(2\omega \tau \right),\\ \epsilon \omega =\epsilon k\sin \left(\omega \tau \right)+\epsilon k\sin \left(2\omega \tau \right),\end{array}$ (2.6)

By computation, we get

$\cos \left(\omega \tau \right)=-\frac{2{\omega }^2-3}{4{\epsilon }^2{\omega }^2+4{\omega }^4-8{\omega }^2+3},\sin \left(\omega \tau \right)=\frac{2\epsilon \omega }{4{\epsilon }^2{\omega }^2+4{\omega }^4-8{\omega }^2+3}$ (2.7)

Solve $\omega$ from Equations (2.6) to get

$4{\epsilon }^4{\omega }^2+8{\epsilon }^2{\omega }^4-16{\epsilon }^2{\omega }^2+4{\omega }^6+5{\epsilon }^2-16{\omega }^4+21{\omega }^2-9=0$ (2.8)

Therefore, Hopf line is given by the curves

$\tau =\frac{1}{\omega }\arccos \left(\frac{-\left(2{\omega }^2-3\right)}{4{\epsilon }^2{\omega }^2+4{\omega }^4-8{\omega }^2+3}\right)$

and Triple-zero bifurcation happens at $\left(\epsilon ,k,\tau \right)=\left(1.3410,0.3729,0.8946\right)$ given that (H₁): The characteristic Equation (2.5) has no eigenvalues $\lambda =i\omega \left(\omega >0\right)$ ;

(H₂): The characteristic Equation (2.5) has zero eigenvalue with algebra multiplicity 3 and geometry multiplicity 1.

3. Triple-Zero Singularity in System

By parameter perturbation, Equation (1.1) is written as

$\begin{array}{l}\ddot{x}+{\epsilon }_0\left(x^2-1\right)\stackrel{\dot{}}{x}+x\\ ={\epsilon }_0{k}_{0}x\left(t-\tau \right)+{\epsilon }_0{k}_{0}x\left(t-2\tau \right)+\left(\epsilon -{\epsilon }_0\right)\left(k_{0}x\left(t-\tau \right)+k_{0}x\left(t-2\tau \right)\right)\\ +{\epsilon }_0\left(k-k_0\right)\left(x\left(t-\tau \right)+x\left(t-2\tau \right)\right)-\left(\epsilon -{\epsilon }_0\right)\left(x^2-1\right)\stackrel{\dot{}}{x}-\beta {\stackrel{\dot{}}{x}}^2+\gamma x^4,\end{array}$ (3.1)

for any $x\in R$ . By Rieze representation theorem, the linearized version of system (3.1) can be written as

$\stackrel{\dot{}}{X}={\displaystyle {\int }_{-2\tau }^0}\text{d}\eta \left(\theta \right)X_t\left(\theta \right)$ (3.2)

with $X={\left(x\mathrm{,}\stackrel{\dot{}}{x}\right)}^{\text{T}}$ . For any $\varphi \in C\left(\left[-2\tau \mathrm{,0}\right]\mathrm{,}{R}^2\right)$ , define the linear operator L by

$L\varphi =A\varphi \left(0\right)+B\varphi \left(-\tau \right)+C\varphi \left(-2\tau \right),\theta =0$ (3.3)

with

$A=\left(\begin{array}{cc}0& 1\\ -1& {\epsilon }_0\end{array}\right),B=\left(\begin{array}{cc}0& 0\\ {\epsilon }_0{k}_0& 0\end{array}\right),C=\left(\begin{array}{cc}0& 0\\ {\epsilon }_0{k}_0& 0\end{array}\right).$ (3.4)

and introduce the linear operator $L_v$ via small parameter perturbation to get

$L_v=A_v\varphi \left(0\right)+B_v\varphi \left(-\tau \right)+C_v\varphi \left(-2\tau \right)$ (3.5)

with

$\begin{array}{l}{A}_v=\left(\begin{array}{cc}0& 0\\ 0& \epsilon -{\epsilon }_0\end{array}\right),B_v=\left(\begin{array}{cc}0& 0\\ \left(\epsilon -{\epsilon }_0\right)k_0+{\epsilon }_0\left(k-k_0\right)& 0\end{array}\right)\\ C_v=\left(\begin{array}{cc}0& 0\\ \left(\epsilon -{\epsilon }_0\right)k_0+{\epsilon }_0\left(k-k_0\right)& 0\end{array}\right)\end{array}$ (3.6)

Based on the foundation theory of DDEs, the linear operator L defines a strong continuous semigroup with infinitesimal generator $\mathcal{A}$ to be expressed as

$\mathcal{A}\varphi =\{\begin{array}{l}\frac{\text{d}\varphi }{\text{d}\theta },-2\tau \le \theta <0,\hfill \\ L\varphi \left(\theta \right),\theta =0\hfill \end{array}$ (3.7)

The adjoint operator $A^{\mathrm{<mo>\; *\; </mo>}}$ in its conjugate phase space $C\left(\left[\mathrm{0,2}\tau \right]\to R^2\right)$ is written as

${\mathcal{A}}^{*}\psi =\{\begin{array}{l}-\frac{\text{d}\psi }{\text{d}s},0<s\le 2\tau ,\hfill \\ L^{*}\psi \left(s\right),s=0\hfill \end{array}$ (3.8)

with

$L^{*}\psi \left(s\right)={\displaystyle {\int }_{-2\tau }^0}\text{d}{\eta }^{\text{T}}\left(-s\right)\psi \left(s\right)=A^{\text{T}}\psi \left(0\right)+B^{\text{T}}\psi \left(\tau \right)+C^{\text{T}}\psi \left(2\tau \right)$ (3.9)

For any $\varphi \in C\mathrm{,}\psi \in C^{\mathrm{<mo>\; *\; </mo>}}$ , define the inner product as the following bilinear form

$\langle \psi \mathrm{,}\varphi \rangle ={\psi }^{\text{T}}\left(0\right)\varphi \left(0\right)-{\displaystyle {\int }_{-2\tau }^0}{\displaystyle {\int }_0^{\theta }}{\psi }^{\text{T}}\left(\xi -\theta \right)\text{d}\eta \left(\theta \right)\varphi \left(\xi \right)\text{d}\xi$ (3.10)

Bogokonav and Takens firstly put forward that Hopf and homoclinic bifurcating from BT point in ODEs. As for DDEs, authors have also verified nongenerate TB bifurcation branches and homoclinic connection while triple zero bifurcation occurs. The following theorem defines triple-zero singularity.

Theorem 1 [6] Under assumption (H₁) and (H₂), the linearized version of DDEs (3.1) has T-B singularity if and only if the following condition holds:

1) $rank\left(A+B+C\right)=n-1$ ;

2) if $N\left(A+B+C\right)=span\left\{q_0\right\}$ , then $q_0\in R^2$ and $\left(B\tau +2C\tau +I\right)q_0\in Range\left(A+B+C\right)$ ;

3) if $\left(A+B+C\right)q_1=\left(B\tau +I+2C\tau \right)q_0$ , then $q_1\in R^2$ and $\left(B\tau +2C\tau +I\right)q_1-\left(\frac{1}{2}B{\tau }^2+2C{\tau }^2\right)q_0\in Range\left(A+B+C\right)$ ;

4) if $\left(A+B+C\right)q_2=\left(B\tau +2C\tau +I\right)q_1-\left(\frac{1}{2}B{\tau }^2+2C{\tau }^2\right)q_0$ , then $q_2\in R^2$ and $\left(\frac{{\tau }^3}{6}B+\frac{8}{6}{\tau }^{3}C\right)q_0-\left(\frac{1}{2}B{\tau }^2+2C{\tau }^2\right)q_1+\left(B\tau +2C\tau +I\right)q_2$ not belongs to $Range\left(A+B+C\right)$ , where ${\varphi }_0\mathrm{,}{\varphi }_2\in R^2$ .

Theorem 1 can be expressed equivalently to

Theorem 2 1) if $N\left(A+B+C\right)=span\left\{q_0\right\}$ , then $rank\left(A+B+C\mathrm{,}\left(B\tau +2C\tau +I\right)q_0\right)=n-1$ ;

2) if $\left(A+B+C\right)q_1=\left(B\tau +I+2C\tau \right)q_0$ , then $rank\left(A+B+C\mathrm{,}\left(B\tau +2C\tau +I\right)q_1-\left(\frac{1}{2}B{\tau }^2+2C{\tau }^2\right)q_0\right)=n-1$ .

3) if $\left(A+B+C\right)q_2=\left(B\tau +2C\tau +I\right)q_1-\left(\frac{1}{2}B{\tau }^2+2C{\tau }^2\right)q_0$ , then $rank\left(A+B+C\mathrm{,}\left(\frac{{\tau }^3}{6}B+\frac{8}{6}{\tau }^{3}C\right)q_0-\left(\frac{1}{2}B{\tau }^2+2C{\tau }^2\right)q_1+\left(B\tau +2C\tau +I\right)q_2\right)=n$ . It is seen that the eigenfunction associated with triple-zero point is chosen as

${\varphi }_0\left(\theta \right)=q_0\mathrm{,}{\varphi }_1\left(\theta \right)=\theta q_0+q_1\mathrm{,}{\varphi }_2\left(\theta \right)=\frac{{\theta }^2}{2}{q}_0+q_1\theta +q_2$

which satisfies

$\mathcal{A}{\varphi }_0\left(\theta \right)=\mathrm{0,}\mathcal{A}{\varphi }_1\left(\theta \right)={\varphi }_0\left(\theta \right),\mathcal{A}{\varphi }_2\left(\theta \right)={\varphi }_1(\; \theta \; )$

since $q_0=\left(\begin{array}{c}1\\ 0\end{array}\right)$ and $q_1=\left(\begin{array}{c}0\\ 1\end{array}\right)$ , $q_2=\left(\begin{array}{c}1\\ 0\end{array}\right)$ satisfies Theorem 1 and Theorem 2, the triple-zero singularity is verified.

The eigenfunction of the adjoint operator $A^{\mathrm{<mo>\; *\; </mo>}}$ is defined and given that

${\mathcal{A}}^{*}{p}_2\left(s\right)=0,{\mathcal{A}}^{*}{p}_1\left(s\right)=p_2\left(s\right),{\mathcal{A}}^{*}{p}_0\left(s\right)=p_1(\; s\; )$

with $p_0\mathrm{,}{p}_1\mathrm{,}{p}_2\in C^{\mathrm{<mo>\; *\; </mo>}}$ , and we also have the following theorem

Theorem 3 The eigenfunction of the adjoint operator $A^{\mathrm{<mo>\; *\; </mo>}}$ can be expressed as $p_2\left(s\right)={\psi }_2$ and $p_1\left(s\right)=-s{\psi }_2+{\psi }_1$ , $p_0\left(s\right)=\frac{s^2}{2}{\psi }_2-s{\psi }_1+{\psi }_0$ , with $0\le s\le 2\tau$ , which satisfies

1) ${\psi }_2^{\text{T}}\left(A+B+C\right)=0$ ;

2) ${\psi }_1^{\text{T}}\left(A+B+C\right)={\psi }_2^{\text{T}}\left(I+\tau B+2\tau C\right)$ ;

3) ${\psi }_1^{\text{T}}\left(I+\tau B+2\tau C\right)-{\psi }_2^{\text{T}}\left(\frac{{\tau }^2}{2}B+2{\tau }^{2}C\right)={\psi }_0^{\text{T}}\left(A+B+C\right)$ ;

4)

${\psi }_2^{\text{T}}{q}_2+{\psi }_{2}B{q}_0\frac{{\tau }^3}{6}-{\psi }_{2}B{q}_1\frac{{\tau }^2}{2}+{\psi }_{2}B{q}_2\tau +{\psi }_{2}C{q}_0\frac{8{\tau }^3}{6}-{\psi }_{2}C{q}_{1}2{\tau }^2+{\psi }_{2}C{q}_{2}2\tau =1$ ;

5) $\begin{array}{l}{\psi }_1^{\text{T}}ph{i}_1+\frac{{\tau }^3}{6}{\psi }_1^{\text{T}}B{q}_0-\frac{{\tau }^2}{2}{\psi }_2^{\text{T}}B{q}_1-\frac{{\tau }^2}{2}{\psi }_1^{\text{T}}B{q}_0+\tau {\psi }_1^{\text{T}}B{q}_1+\frac{8{\tau }^3}{6}{\psi }_1^{\text{T}}C{q}_0\\ -2{\tau }^2{\psi }_2^{\text{T}}C{q}_1-2{\tau }^2{\psi }_1^{\text{T}}C{q}_0+2\tau {\psi }_1^{\text{T}}C{q}_1=1\end{array}$ ;

6)

${\psi }_0^{\text{T}}{q}_0+\frac{{\tau }^3}{6}{\psi }_{2}B{q}_0-\frac{{\tau }^2}{2}{\psi }_{1}B{q}_0+\tau {\psi }_{0}B{q}_0+\frac{8{\tau }^3}{6}{\psi }_{2}C{q}_0-2{\tau }^2{\psi }_{1}C{q}_0+2\tau {\psi }_{0}C{q}_0=1$ ;

7) ${\psi }_2{q}_1-{\psi }_{2}B{q}_0\frac{{\tau }^2}{2}+{\psi }_{2}B{q}_1\tau -{\psi }_{2}C{q}_{0}2{\tau }^2+{\psi }_{2}C{q}_{1}2\tau =0$ .

Choose ${\psi }_0=\left(\begin{array}{c}-2.0292\\ 0.8668\end{array}\right)$ and ${\psi }_1=\left(\begin{array}{c}-4.3627\\ 1.8636\end{array}\right)$ , ${\psi }_2=\left(\begin{array}{c}-2.4932\\ 1.8592\end{array}\right)$ then it is verified that conditions (1) - (7) are satisfied in theorem 3.

4. Triple-Zero Bifurcation of Reduction on Center Manifold

Define P is projection operator from phase space to its eigenspace associated with triple zero singularity. Assume $\Lambda =\left\{0,0\right\}$ and $P_{\Lambda }$ is the related invariant eigenspace and $Q_{\Lambda }$ is the complementary subspace of $P_{\Lambda }$ . Now we consider the decomposition of the phase space as $C=P_{\Lambda }\oplus Q_{\Lambda }$ . Define the dual eigenspace of $P_{\Lambda }$ as $P_{\Lambda }^{\mathrm{<mo>\; *\; </mo>}}=span\left\{p_0\left(s\right)\mathrm{,}{p}_1\left(s\right)\mathrm{,}{p}_2\left(s\right)\right\}$ , then it is directly derived that

${\mathcal{A}}^{*}{p}_2\left(s\right)=0,{\mathcal{A}}^{*}{p}_1\left(s\right)=p_2\left(s\right),{\mathcal{A}}^{*}{p}_0\left(s\right)=p_1(\; s\; )$

To apply the method of Faria and Magalhães [15] [16] , it is necessary to enlarge the phase space C to its extended phase space BC which is uniformly continuous on $\left[-2\tau \mathrm{,0}\right]$ however may have a jump discontinuity at 0. As referred above to have the decomposition of phase space as $C=P_{\Lambda }\oplus Q$ , the projection $\Pi \mathrm{<mo>\; :\; </mo>}BC\to P_{\Lambda }$ leads to the decomposition $BC=P_{\Lambda }\oplus Ker\left(\Pi \right)$ . Now, by setting $X_t=y_1{\varphi }_0+y_2{\varphi }_1+y_3{\varphi }_2+v_t$ , with $v_t\in Ker\left(\Pi \right)$ , with the fundamental theory of DDEs, the projection of DDEs (3.1) is written as

#Math_98# (4.2)

with $J=\left(\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 0& 0& 0\end{array}\right),\Psi \left(0\right)=col\left({\psi }_0,{\psi }_1,{\psi }_2\right)$ , and also the bifurcation equation

${v^{\prime }}_t=\mathcal{A}{v}_t+\left(I-\Pi \right)X_{0}F\left(y_1{\varphi }_0+y_2{\varphi }_1+y_3{\varphi }_2+v_t\right)\mathrm{,}$ (4.3)

with the definition

$X_0=\{\begin{array}{l}I\theta =0\hfill \\ 0-2\tau \le \theta \le 0\hfill \end{array}$ (4.4)

herein

$\begin{array}{l}F\left(y_1{\varphi }_0+y_2{\varphi }_1+y_3{\varphi }_2+v_t\right)\\ =\left(\begin{array}{c}0\\ -{\epsilon }_0{\left(y_1+y_3+v_t^{\left(1\right)}\left(0\right)\right)}^2\left(y_2+v_t^{\left(2\right)}\left(0\right)\right)-\beta {\left(y_2+v_t^{\left(2\right)}\left(0\right)\right)}^2+\gamma {\left(y_1+y_3+v_t^{\left(1\right)}\left(0\right)\right)}^4\end{array}\right)\end{array}$ (4.5)

Set

$v_t=\left(\begin{array}{c}{\displaystyle {\sum }_{i+j+k=2}{g}_{ijk}^{\left(1\right)}{y}_1^i{y}_2^j{y}_3^k}+\cdots \\ {\displaystyle {\sum }_{i+j+k=2}{g}_{ijk}^{\left(2\right)}{y}_1^i{y}_2^j{y}_3^k}+\cdots \end{array}\right)$ (4.6)

Accordingly, setting $\epsilon -{\epsilon }_0=\mu v,k-k_0=sv,\tau -{\tau }_0=mv$ , the Taylor expansion of Equation (4.2) is rewritten as sum of series as

$\begin{array}{l}{\Psi }^{\text{T}}\left(0\right)\left(L_v\left(y_1{\varphi }_0+y_2{\varphi }_1+y_3{\varphi }_2+v_t\right)+F\left(y_1{\varphi }_0+y_2{\varphi }_1+y_3{\varphi }_2+v_t\right)\right)\\ ={\displaystyle {\sum }_{i+j+k+s\ge 2}{f}_{ijks}^{\left(1\right)}{y}_1^i{y}_2^j{y}_3^k{v}^s}\end{array}$ (4.7)

and

$\left(I-\Pi \right)X_{0}F\left(y_1{\varphi }_0+y_2{\varphi }_1+y_3{\varphi }_2+v_t\right)={\displaystyle {\sum }_{i+j+k+s\ge 2}{f}_{ijks}^{\left(2\right)}{y}_1^i{y}_2^j{y}_3^k{v}^s}$ (4.8)

Let $H_j^4$ denote the homogeneous polynomial space of degree $j\left(j\ge 2\right)$ with respect to variable $\left(y_1\mathrm{,}{y}_2\mathrm{,}{y}_3\mathrm{,}v\right)$ . To do further calculation of change of variable $\left(y_1\mathrm{,}{y}_2\mathrm{,}{y}_3\mathrm{,}v\right)$ , one defines the conjugate operator $M_j^1$ in $H_j^4$ as

$M_j^{1}p\left(y\mathrm{,}v\right)=p_y\left(y\mathrm{,}v\right)\left(Jy\right)-Jp\left(y\mathrm{,}v\right)$ (4.9)

On the center manifold, the normal form of Equations (4.5) is expressed by

$\left(\begin{array}{c}{y^{\prime }}_1\\ {y^{\prime }}_2\\ {y^{\prime }}_3\end{array}\right)=B\left(\begin{array}{c}{y}_1\\ y_2\\ y_3\end{array}\right)+g_2^{\left(1\right)}\left(y_1\mathrm{,}{y}_2\mathrm{,}{y}_3\mathrm{,0}\right)+g_3^{\left(1\right)}\left(y_1\mathrm{,}{y}_2\mathrm{,}{y}_3\mathrm{,0}\right)$

with

$\begin{array}{l}{g}_2^{\left(1\right)}\left(y_1\mathrm{,}{y}_2\mathrm{,}{y}_3\mathrm{,0}\right)=Pro{j}_{Im{\left(M_2^{\left(1\right)}\right)}^c}{f}_2^{\left(1\right)}\left(y_1\mathrm{,}{y}_2\mathrm{,}{y}_3\mathrm{,0}\right)\\ g_3^{\left(1\right)}\left(y_1\mathrm{,}{y}_2\mathrm{,}{y}_3\mathrm{,0}\right)=Pro{j}_{Im{\left(M_3^{\left(1\right)}\right)}^c}{\tilde{f}}_3^{\left(1\right)}\left(y_1\mathrm{,}{y}_2\mathrm{,}{y}_3\mathrm{,0}\right)\end{array}$

The bases of $Im\left(M_2^1\right)$ is expressed as

$\begin{array}{l}\left(\begin{array}{c}{y}_1{y}_2\\ 0\\ 0\end{array}\right)\mathrm{,}\left(\begin{array}{c}{y}_2{y}_3\\ 0\\ 0\end{array}\right)\mathrm{,}\left(\begin{array}{c}{y}_1{y}_3+y_2^2\\ 0\\ 0\end{array}\right)\mathrm{,}\left(\begin{array}{c}{y}_3^2\\ 0\\ 0\end{array}\right)\mathrm{,}\left(\begin{array}{c}{y}_{2}v\\ 0\\ 0\end{array}\right)\mathrm{,}\\ \left(\begin{array}{c}{y}_{3}v\\ 0\\ 0\end{array}\right)\mathrm{,}\left(\begin{array}{c}-y_1^2\\ 2y_1{y}_2\\ 0\end{array}\right),\left(\begin{array}{c}-y_2^2\\ 2y_2{y}_3\\ 0\end{array}\right)\mathrm{,}\left(\begin{array}{c}-y_1{y}_2\\ y_1{y}_3+y_2^2\\ 0\end{array}\right)\mathrm{,}\left(\begin{array}{c}-y_1{y}_3\\ y_2{y}_3\\ 0\end{array}\right)\mathrm{,}\\ \left(\begin{array}{c}-y_2{y}_3\\ y_3^2\\ 0\end{array}\right)\mathrm{,}\left(\begin{array}{c}-y_{1}v\\ y_{2}v\\ 0\end{array}\right)\mathrm{,}\left(\begin{array}{c}-y_{2}v\\ y_{3}v\\ 0\end{array}\right)\mathrm{,}\left(\begin{array}{c}0\\ -y_1^2\\ 2y_1{y}_2\end{array}\right)\mathrm{,}\left(\begin{array}{c}0\\ -y_2^2\\ 2y_2{y}_3\end{array}\right)\mathrm{,}\\ \left(\begin{array}{c}0\\ y_3^2\\ 0\end{array}\right)\mathrm{,}\left(\begin{array}{c}0\\ -y_1{y}_2\\ y_1{y}_3+y_2^2\end{array}\right)\mathrm{,}\left(\begin{array}{c}0\\ -y_2{y}_3\\ y_3^2\end{array}\right)\mathrm{,}\left(\begin{array}{c}0\\ -y_1{y}_3\\ y_2{y}_3\end{array}\right)\mathrm{,}\left(\begin{array}{c}0\\ -y_{1}v\\ y_{2}v\end{array}\right)\mathrm{,}\\ \left(\begin{array}{c}0\\ -y_{2}v\\ y_{3}v\end{array}\right)\mathrm{,}\left(\begin{array}{c}0\\ y_{3}v\\ 0\end{array}\right)\end{array}$

The bases of the $Ker\left(M\right)$ is listed as

$\left(\begin{array}{c}{y}_1^2\\ y_1{y}_2\\ \frac{y_2^2}{2}\end{array}\right)\mathrm{,}\left(\begin{array}{c}{y}_1^2\\ y_1{y}_2\\ y_1{y}_3\end{array}\right)\mathrm{,}\left(\begin{array}{c}0\\ y_1^2\\ y_1{y}_2\end{array}\right)\mathrm{,}\left(\begin{array}{c}0\\ 0\\ y_1^2\end{array}\right)\mathrm{,}\left(\begin{array}{c}{y}_{1}v\\ y_{2}v\\ y_{3}v\end{array}\right)\mathrm{,}\left(\begin{array}{c}0\\ y_{1}v\\ y_{2}v\end{array}\right)\mathrm{,}\left(\begin{array}{c}0\\ 0\\ y_{1}v\end{array}\right)$

Based on the projection operator $M_j^{1}p\left(y\mathrm{,0}\right)$ with decomposition of homogeneous polynomial space $H^j$ , Equations (3.1) is transformed into its simplified version

$\begin{array}{l}{y^{\prime }}_1=y_2,\\ {y^{\prime }}_2=y_3,\\ {y^{\prime }}_3=a_0{y}_{3}v+a_1{y}_{2}v+a_2{y}_{1}v+a_3{y}_1{y}_3+\frac{a_4}{2}{y}_2^{(2)}+a_5{y}_1{y}_2+\frac{a_6}{4}{y}_1^2,\end{array}$ (4.10)

with

$\begin{array}{l}{a}_0={\psi }_2^{\left(2\right)}\mu +2{\psi }_2^{\left(2\right)}{k}_0\mu +2{\psi }_2^{\left(2\right)}s{\epsilon }_0+\frac{5}{2}{\psi }_2^{\left(2\right)}{k}_0\mu {\tau }_0^2+\frac{5}{2}{\psi }_2^{\left(2\right)}s{\tau }_0^2{\epsilon }_0\\ +8{\psi }_2^{\left(2\right)}\mu vs+10{\psi }_2^{\left(2\right)}\mu vs{\tau }_0^2+5{\psi }_2^{\left(2\right)}{k}_{0}m{\tau }_0{\epsilon }_0-12\mu {\psi }_1^{\left(2\right)}s{\tau }_{0}v-3k_{0}m{\psi }_1^{\left(2\right)}{\epsilon }_0\\ -3k_0\mu {\psi }_1^{\left(2\right)}{\tau }_0-3{\psi }_1^{\left(2\right)}s{\tau }_0{\epsilon }_0+8\mu {\psi }_0^{\left(2\right)}sv+2k_0\mu {\psi }_0^{\left(2\right)}+2{\psi }_0^{\left(2\right)}s{\epsilon }_0+\mu {\psi }_0^{\left(2\right)},\\ a_1=-12\mu {\psi }_2^{\left(2\right)}s{\tau }_{0}v-3k_{0}m{\psi }_2^{\left(2\right)}{\epsilon }_0-3k_0\mu {\psi }_2^{\left(2\right)}{\tau }_0+8\mu {\psi }_1^{\left(2\right)}sv\\ -3{\psi }_2^{\left(2\right)}s{\tau }_0{\epsilon }_0+2k_0\mu {\psi }_1^{\left(2\right)}+2{\psi }_1^{\left(2\right)}s{\epsilon }_0+\mu {\psi }_1^{\left(2\right)},\end{array}$

$\begin{array}{l}{a}_2=8\mu {\psi }_2^{\left(2\right)}sv+2k_0\mu {\psi }_2^{\left(2\right)}+2{\psi }_2^{\left(2\right)}s{\epsilon }_0+\mu {\psi }_2^{\left(2\right)},\\ a_3=-2\beta {\psi }_0^{\left(2\right)}-2\beta {\psi }_2^{\left(2\right)},\\ a_4=-2{\psi }_0^{\left(2\right)}\beta ,\\ a_5=-2{\psi }_1^{\left(2\right)}\beta ,\\ a_6=2{\psi }_2^{\left(2\right)}\beta \end{array}$

By Equation (4.10), one derives that both double-zero bifurcation line and Hopf-bifurcation line are produced as varying variable v and two lines are merged at $v=0$ . Accordingly, the produced fold line is expressed as:

$\left\{\left(\mu ,s,m\right)|8\mu {\psi }_2^{\left(2\right)}s{v}^2+2k_0\mu {\psi }_2^{\left(2\right)}v+2{\psi }_2^{\left(2\right)}sv{\epsilon }_0+\mu {\psi }_2^{\left(2\right)}v=0\right\}$ (4.11)

The double-zero bifurcation line is denoted as

$\begin{array}{l}\{\left(\mu \mathrm{,}s\mathrm{,}m\right)\mathrm{<mo>\; |\; </mo>}-12\mu {\psi }_2^{\left(2\right)}s{\tau }_0{v}^2-3k_{0}m{\psi }_2^{\left(2\right)}v{\epsilon }_0-3k_0\mu {\psi }_2^{\left(2\right)}{\tau }_{0}v+8\mu {\psi }_1^{\left(2\right)}s{v}^2\\ -3{\psi }_2^{\left(2\right)}s{\tau }_{0}v{\epsilon }_0+2k_0\mu {\psi }_1^{\left(2\right)}v+2{\psi }_1^{\left(2\right)}sv{\epsilon }_0+\mu {\psi }_1^{\left(2\right)}v=0\}\end{array}$ (4.12)

and Hopf line is also expressed as

$\begin{array}{l}{\omega }^2=12\mu {\psi }_2^{\left(2\right)}s{\tau }_0{v}^2+3k_{0}m{\psi }_2^{\left(2\right)}v{\epsilon }_0+3k_0\mu {\psi }_2^{\left(2\right)}{\tau }_{0}v-8\mu {\psi }_1^{\left(2\right)}s{v}^2+3{\psi }_2^{\left(2\right)}s{\tau }_{0}v{\epsilon }_0\\ -2k_0\mu {\psi }_1^{\left(2\right)}v-2{\psi }_1^{\left(2\right)}sv{\epsilon }_0-\mu {\psi }_1^{\left(2\right)}v,\end{array}$

$\begin{array}{l}\{\left(\mu ,s,m\right)|-{\omega }^{2}v{\psi }_2^{\left(2\right)}\mu -{\omega }^{2}v\mu {\psi }_0^{\left(2\right)}+2v{\psi }_2^{\left(2\right)}{k}_0\mu +2v{\psi }_2^{\left(2\right)}s{\epsilon }_0+8{\psi }_2^{\left(2\right)}\mu v^{2}s\\ -\frac{5}{2}{\omega }^{2}v{\psi }_2^{\left(2\right)}{k}_0\mu {\tau }_0^2-\frac{5}{2}{\omega }^{2}v{\psi }_2^{\left(2\right)}s{\tau }_0^2{\epsilon }_0-10{\omega }^2{\psi }_2^{\left(2\right)}\mu v^{2}s{\tau }_0^2+12{\omega }^2\mu {\psi }_1^{\left(2\right)}s{\tau }_0{v}^2\\ +3{\omega }^{2}v{k}_{0}m{\psi }_1^{\left(2\right)}{\epsilon }_0+3{\omega }^{2}v{k}_0\mu {\psi }_1^{\left(2\right)}{\tau }_0+3{\omega }^{2}v{\psi }_1^{\left(2\right)}s{\tau }_0{\epsilon }_0-5{\omega }^{2}v{\psi }_2^{\left(2\right)}{k}_{0}m{\tau }_0{\epsilon }_0\\ -2{\omega }^{2}v{\psi }_2^{\left(2\right)}{k}_0\mu -2{\omega }^{2}v{\psi }_2^{\left(2\right)}s{\epsilon }_0-8{\omega }^2{\psi }_2^{\left(2\right)}\mu v^{2}s-8{\omega }^2\mu {\psi }_0^{\left(2\right)}s{v}^2-2{\omega }^{2}v{k}_0\mu {\psi }_0^{\left(2\right)}\\ -2{\omega }^{2}v{\psi }_0^{\left(2\right)}s{\epsilon }_0+v{\psi }_2^{\left(2\right)}\mu =0\}\end{array}$ (4.13)

Set $m=0$ , then double zero line which is drawn in blue color is in coincidence with the saddle node line expressed by the above parameter sets. The bifurcation of triple zero singularity is acquired as Hopf lines merged with double zero lines, as shown in Figure 1(a). Correspondingly, both Hopf line and double zero bifurcation line of van der Pol system (1.1) are drawn in Figure 1(b), wherein eigenvalues of zero real part are determined by Equation (2.6). In addition, by the application of DDE-Biftool software, the bifurcating of homoclinic line of system (4.10) and van der Pol system (1.1) are also respectively produced as varying parameter pairs $\left(\epsilon \mathrm{,}k\right)$ continuously at fixed time delay $\tau =0.8946$ , as shown in Figure 2(a) and Figure 2(b). At fixed time delay $\tau =1.094$ , the complex dynamics near triple zero point is also simulated, as shown in Figures 3(a)-(f). The scenarios of period doubling bifurcation of limit cycle are exhibited

Figure 1. The bifurcation of triple zero singularity. (a) Bifurcation lines of system (4.10) determined by parameter bifurcation sets with $m=0$ ; (b) The associated bifurcation lines of van der Pol system (1.1) with fixed value $\tau =0.8946$ .

Figure 2. The homoclinic loops produced near triple zero bifurcation point. (a) The homoclinic solution of Equation (4.10) for $\mu =-0.0643,s=0.0329$ ; (b) The continuous homoclinic solutions of system (1.1).

Figure 3. The phase portrait of the bifurcating limit cycles near the triple zero bifurcation point of van der Pol system (1.1) with $\tau =1.094$ . (a) The double period limit cycle; (b) Period-3 limit cycle; (c) Period-4 limit cycle; (d) Period-8 limit cycle; (e) quasi-period solution; (f) Chaos attractor.

Figure 4. The scenarios of limit cycle bifurcation and the routes to chaos with time delay $\tau =1.094$ . (a) The bifurcation on Poincare section with chosen $\epsilon =0.7593$ ; (b) The bifurcation on Poincare section with chosen $k=0.475$ .

on phase space as either k descending or $\epsilon$ ascending, as shown in Figure 4(a) and Figure 4(b) respectively, and the routes to chaos are also simulated by chosen Poincare section $x\left(t-\tau \right)=0$ . The period doubling bifurcation brings forth periodic solutions oscillating two times within one prolonging period instantaneously [19] [20] .

5. Conclusion

As the double zero bifurcation line merges into Hopf point, the triple zero bifurcation point arose in van der Pol system. The complex dynamics of van der Pol oscillator with delay feedback control was investigated. With the application of analytical technique of dimension reduction method combined with the center manifold theory, the normal form near the triple zero singularity was derived and the bifurcating homoclinic solution was analyzed further. With the aid of technique artificial tool of DDE-Biftool software, the simulated homoclinic solution was plotted in phase space continuously. The complex dynamics as double period limit cycle near bifurcation point of triple zero singularity was observed and the routes to chaos were also simulated by choosing the corresponding Poincare section as varying parameter continuously.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Lin, G. (2007) Periodic Solutions for van der Pol Equation with Time Delay. Applied Mathematics and Computation, 187, 1187-1198. https://doi.org/10.1016/j.amc.2006.09.032
[2]	Jiang, W. and Wei, J. (2008) Bifurcation Analysis in van der Pol’s Oscillator with Delayed Feedback. Journal of Computational and Applied Mathematics, 213, 604-615. https://doi.org/10.1016/j.cam.2007.01.041
[3]	Maccari, A. (2001) The Response of Parametrically Excited van der Pol Oscillator to a Time Delay State Feedback. Nonlinear Dynamics, 26, 105-119.
[4]	Eigoli, A.K. and Khodabakhsh, M. (2011) A Homotopy Analysis Method for Limit Cycle of the van der Pol Oscillator with Delayed Amplitude Limiting, Applied Mathematics and Computation, 217, 9404-9411. https://doi.org/10.1016/j.amc.2011.04.029
[5]	Yu, W. and Cao, J. (2005) Hopf Bifurcation and Stability of Periodic Solutions for van der Pol Equation with Time Delay. Nonlinear Analysis, 62, 141-165. https://doi.org/10.1016/j.na.2005.03.017
[6]	Xu, Y. and Huang, M. (2008) Homoclinic Orbits and Hopf Bifurcations in Delay Differential Systems with T-B Singularities. Journal of Differential Equations, 244, 582-598. https://doi.org/10.1016/j.jde.2007.09.003
[7]	Liao, X. (2005) Hopf and Resonant Codimension Two Bifurcation in van der Pol Equation with Two Time Delays. Chaos Solitons Fractals, 23, 857-871. https://doi.org/10.1016/j.chaos.2004.05.048
[8]	Liu, X. and Wang, J.L. (2013) Bogdanov-Takens and Triple Zero Bifurcations of a Delayed Modified Leslie-Gower Predator Prey System. Abstract and Applied Analysis, 2013, Article ID: 605471. https://doi.org/10.1155/2013/605471
[9]	Campell, S.A. and Yuan, Y. (2008) Zero Sigularities of Codimension Two and Three in Delay Differential Euqations. Nonlinearity, 21, 2671-2691. https://doi.org/10.1088/0951-7715/21/11/010
[10]	He, X., Li, C. and Shu, Y.L. (2012) Triple-Zero Bifurcation in van der Pol’s Oscillator with Delayed Feedback. Communications in Nonlinear Science and Numerical Simulation, 17, 5229-5239. https://doi.org/10.1016/j.cnsns.2012.05.001
[11]	Li, Y. and Jiang, W. (2011) Hopf and Bogdanov-Takens Bifurcation in a Coupled Fitzhugh-Nagumo Neural System with Time Delay. Nonlinear Dynamics, 65, 161-173. https://doi.org/10.1007/s11071-010-9881-5
[12]	Hoff, A., dos Santos, J.V., et al. (2004) Numerical Bifurcation Analysis of Two Coupled Fitahugh-Nagumo Oscillators. The European Physical Journal B, 87, Article No. 151. https://doi.org/10.1140/epjb/e2014-50170-9
[13]	Kuznetsov, Y.A. (2013) Elements of Applied Bifurcation Theory. Springer Science & Business Media, Berlin.
[14]	Faria, T. and Magalhaes, L.T. (1995) Normal Forms for Retarded Functional Differential Equations and Applications to Bogdanov-Takens Sigularity. Journal of Differential Equations, 12, 201-224. https://doi.org/10.1006/jdeq.1995.1145
[15]	Faria, T. and Magalhaes, L.T. (1995) Normal Forms for Retarded Functional Equations with Parameters and Applications to Hopf Bifurcation. Journal of Differential Equations, 122, 181-200. https://doi.org/10.1006/jdeq.1995.1144
[16]	Guckenheimer, J. and Holmes, P. (1983) Nonlinear Oscillations, Dynamical Systems, and Bifurcation of Vector Fields. Springer, New York. https://doi.org/10.1007/978-1-4612-1140-2
[17]	Engelborghs, K., Luzyanina, T. and Roose, D. (2002) Numerical Bifurcation Analysis of Delay Differentia, Equations Using DDE-BIFTOOL. ACM Transactions on Mathematical Software (TOMS), 28, 1-21. https://doi.org/10.1145/513001.513002
[18]	Lin, Y. (2014) Stability and Oscillation of the Solutions for a Coupled FHN Model with Time Delays. In: Intelligent Computing in Bioinformatics, Springer, Berlin, 164-174. https://doi.org/10.1007/978-3-319-09330-7_21
[19]	Sun, J.Q. and Ding, G. (2013) Advances in Analysis and Control of Time-Delayed Dynamical Systems. World Scientific, Singapore. https://doi.org/10.1142/8878
[20]	Flunkert, V., Fisher, I. and Scröll, E. (2013) Dynamics, Control and Information in Delayed-Coupled Systems. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 371, 1-4. https://doi.org/10.1098/rsta.2012.0465