Complex dynamics and Bogdanov-Takens bifurcations in a retarded van der Pol-Duffing oscillator with positional delayed feedback

1.   Introduction

It is well known that the dynamics of harmonic oscillators has long been studied for numerous applications in electrodynamics, engineering, electronics, neurology and biological systems. Time-delayed feedback control has been used to a variety of disciplines of inquiry, including chaos communication, electronic systems, and engineering. It is known that the amplitude of oscillations always differs for the positional delayed feedback. The purpose of the present study is to investigate the global dynamics of a retarded van der-Pol Duffing oscillator with retarded damping and delayed position feedback, which is a new idea in the study of harmonic oscillators. Due to the numerous uses of the harmonic oscillator, including in physics and many other disciplines, it has been discussed for a very long period. If the frictional force is proportional to the velocity, the harmonic oscillator is referred to as a damped oscillator. In the case of damped oscillators, the vibration amplitude decreases with time. Damping is important in realistic oscillatory systems, which can be positive or negative.

We get delay differential equation by introducing time delay into the ordinary differential equation. Clearly it is more realistic. Following ground-breaking discovery done by Pyragas ^[1], time-delayed feedback control has been applied to numerous areas of study, including chaos communication, engineering and electronic systems. From the study of Atay ^[2], we know that positional delayed feedback change the amplitude of oscillations. We can describe the damped harmonic oscillator with delayed feedback by the following functional differential equation:

where $x(t)\in\mathbf{R}$ denote the displacement from the point of oscillation, $a > 0$ denote the stiffness of the spring element, $g$ is a $C^r\; (r \geq3)$ smooth function which narrates the delayed feedback, $\epsilon$ is the damping coefficient, $\eta$ denote the time delay, and $b$ is the Duffing coefficient of order three. Equation (1.1) has been analyzed by many authors (cf. Jiang and Song ^[3], Campbell et al. ^[4], Song et al. ^[5], Song and Xu ^[6], Cao et al. ^[7], Campbell and Yuan ^[8]).

So far, there have been numerous research findings on Hopf bifurcation, Hopf-zero bifurcation and triple-zero bifurcation of several oscillators studied by Cao and Yuan ^[9], Qiao et al. ^[10]. But, there are limited number of research articles that described the B-T bifurcation of the harmonic oscillators with multiple delay and retarded damping. Inspired from the discussions in ^{[11,12,13,14]} time delay is incorporated in Eq (1.1), which then extended to the following retarded van der Pol-Duffing equation with positional delayed feedback:

It is simple to verify that the system's origin Eq (1.2) is a double-zero singularity subjected to some parametric conditions. The primary target of this study is to discuss the B-T bifurcation and extract the related normal form of the dynamical Eq (1.2). We know that B-T singularity is a critical point at which the Jacobian matrix has double zero eigenvalue and geometric multiplicity one. Also, we know that the order of the ordinary differential equation is at least two if it contains a B-T singularity, which is not true for the delay differential equation.

The system of delay differential equations can be converted to a two-dimensional ordinary differential equation centred on the B-T singular point, using the center manifold reduction and normal form theory. In fact, a variety of dynamical systems, including chemical reaction and neural network models, can accurately depict the B-T bifurcation under certain critical circumstances. This paper will investigate the B-T bifurcation of a retarded delay differential equation with damping effect.

The remaining part of this article is outlined as follows: The parametric restrictions under which the origin is a B-T singularity are stated in Section 2. Section 3 provides a detailed description of the B-T singularity as well as the second and third order normal forms corresponding bifurcation. In Section 4, we carried out numerical simulations to verify the analytical findings made in this study. In Section 5, the article concludes with a brief conclusion.

2.   Bogdanov-Takens bifurcation

The Bogdanov-Takens bifurcation is an unusual and fascinating bifurcation of continuous dynamical systems, in which the Jacobian matrix has a double-zero eigenvalue. It's a bit difficult to explain what's going on. For this type of bifurcation, there is a two dimensional parametric space, wherein there exist a curve of saddle-node bifurcations; a curve of Hopf bifurcations; and a curve of homoclinic bifurcations. The intersection of these three curves is know as Bogdanov-Takens bifurcation point. One of the ways that limit cycles are created or annihilated is in the context of a Bogdanov-Takens bifurcation. These global bifurcations are unlike anything we've yet seen in the local bifurcations like Hopf, Saddle-node, Transcritical and Pitchfork bifurcations. Considering the aforementioned discussion, the retarded van der Pol-Duffing oscillator with two delays can represented mathematically with the following retarded delay differential system:

where $\tau$ and $\eta$ represents the retarded delay on damping and positional feedback respectively. Let $x(t) = x_1(t)$ and $\dot{x}_1(t) = x_2(t)$. Then Eq (2.1) becomes

We will study this Eq (2.1) theoretically and numerically for both the focus $(\epsilon > 0)$ and saddle $(b < 0)$ cases.

Let us consider $a > 0, \; g(0) = 0,$ and $g^{\prime}(0) = d.$ Thus, the characteristic equation of the Jacobian matrix of the Eq (2.2) around $(0, 0)$ is given by

Lemma 1. Let $2+\epsilon(\eta-2\tau) > 0, \; d = a$ and $\epsilon = -a\eta$. Then the Eq (2.3) has a double-zero root, namely $\lambda = 0$.

Proof. Here, $E(0) = a-d$. $E(0) = 0$ as $d = a$. Clearly, $E^{'}(\lambda) = 2\lambda+\epsilon e^{-\lambda\tau}-\epsilon\tau\lambda e^{-\lambda \tau}+d\eta e^{-\lambda\eta},$ consequently $E^{'}(0) = 0$ as $\epsilon = -a\eta$. Also, we have $E{''}(0) = 2+\epsilon(\eta-2\tau).$ Thus $E{''}(0)\neq0$, since $2+\epsilon(\eta-2\tau) > 0$. Therefore, we conclude that $0$ is a characteristic root of multiplicity two.

Lemma 2. If $2+\epsilon(\eta-2\tau) > 0, \; d = a$ and $\epsilon = -a\eta$, and $a\in (0, M)$, then all roots of Eq (2.3) except $\lambda = 0$, have non-zero real parts, i.e., origin of the Eq (2.2) is B-T singularity, where $M = \{\min a_j:\; a_j =$ $\frac{\sigma_j^2}{1-\eta \sigma_j\sin{\tau \sigma_j}-\cos{\eta \sigma_j}} >$ $0, \; 1\leq j\leq m, \; and\; \sigma_j$ $\; are\; the\; roots\; of \; the\; equation\; \sigma^2-(\epsilon^2+2a)+$ $\frac{2a\epsilon\sin(\tau-\eta)\sigma}{\sigma} = 0\}.$

Proof. The proof is carried out in the same technique adopted by Sarwardi et al. ^[15].

3.   Normal forms of Bogdanov-Takens bifurcation of the Eq (2.2)

From Lemma 2, we see that the Eq (2.2) experiences B-T bifurcation around the origin, if $\epsilon = -a\eta$ and $d = a\in(0, M)$ hold. Hence, we assume $d$ and $\epsilon$ as bifurcation parameters. Therefore, we may incorporate two small parameters $\nu_1$ and $\nu_2$ in the vicinity of origin and then discuss perturbation of the Eq (2.2). Expanding the Eq (2.2) according to Taylor's theorem, we have

After simplification, on the phase space $C$, we may write the Eq (3.1) as the following retarded functional differential equation :

where $C$ is the Banach space of all continuous functions from $[-\tau, 0]$ to $\mathbf{R}^2$ with supremum norm $|\xi| = \sup_{\theta\in[-\tau, 0]}|\xi(\theta)|,$ $X_t\in C$ is defined by $X_t(\theta) = X(t+\theta)$ for $\theta\in[-\tau, 0],$ and $L(\nu):C\rightarrow \mathbf{R}^2$ is a parameterized family of bounded linear operators defined as follows:

and $H:C\times W\rightarrow \mathbf{R}^2$ is a $C^k\; (k\geq2)$ function satisfying $H (0, 0), \; DH(0, 0) = 0$ with

Now, we will linearize the Eq (3.2) around $(X_t, \nu) = (0, 0)$ and obtained the followings:

where $L(0)(\xi) = A\xi(0)+B_1\xi(-\tau)+B_2\xi(-\eta),$

The linear operator $L_1$ may be defined in the similar fashion.

Now, we will use Riesz representation theorem on the bounded linear operator $L_0$ and obtain a $2\times2$ matrix $\rho(.)$ defined on $[-\tau, 0]$ whose entries are functions of bounded variation, such that

From ^[16] and ^[17], we define the adjoint inner product on $C^*\times C$ as follows:

where $C^* = C([0, \tau], \mathbf{R}^{2*})$ be the adjoint space of $C$. Let $A_0:C\rightarrow C$ be the infinitesimal generator, defined by $A_0(\psi) = \dot{\psi}$. Assuming $\Lambda_0 = \{0, 0\}$, the set of eigenvalues with zero real part of $A_0$, we define the invariant space of $A_0$ associated with $\Lambda_0$, and the dual of it by $P$ and $P^*$ respectively. Now, we will decompose $C$ by the help of adjoint theory for a functional differential equation as follows:

Now, we define the bases for $P$ and $P^*$ by $\Psi$ and $\Phi$ respectively. We are assuming $\Psi$ and $\Phi$ as follows: $\Psi = (\psi_1(\theta), \psi_2(\theta))$ for $-\tau\leq\theta\leq 0$, and $\Phi = (\phi_1(s), \phi_2(s))^T$ for $0\leq s\leq\tau.$ Then $\left < \Phi, \Psi\right > = I_2$, and $\dot{\Psi} = \Psi\bar {B},$ $-\dot{\Phi} = \bar{B}\Phi,$ where $\bar{B} = \left(\begin{array}{cc} 0 & 1 \\ 0 & 0 \\ \end{array} \right).$ After using Lemma $3.1$ of ^[18], we get

and

In the next part of this section, B-T bifurcation of the Eq (3.2) will be discussed. For that we will extend the phase space $C$ by Banach space $BC = \{\psi:[-\tau, 0]\rightarrow \mathbf{R}^2:\psi\; \text{is continuous on}\; [-\tau, 0), \; \text{with a possible jump discontinuity at zero}\}.$ Any element $\psi$ of $BC$ may be assumed as $\psi = \chi+Y_0 e$ and norm of the same is given by $||\psi|| = ||\chi+Y_0c|| = ||\chi||_C+|e|$, where $Y_0$ is a $2\times 2$ matrix valued function defined as follows

Now, in $BC$, the Eq (3.2) may be transformed to the following form:

where $G(w, \nu) = (L(\nu)-L_0)w+H(w, \nu) = L_1(\nu)w+H(w, \nu), \; \nu\in W$ and $\bar{A}:C^1\rightarrow BC$ is the extension of infinitesimal generator $A_0$, it is defined by

Now, we will use the result of ^[19] and letting $w = \mathbf{ \pmb{\mathsf{ ψ}}}(\theta)z(t)+y$ to decompose the Eq (3.10) as follows:

where $\pi:BC\rightarrow P$ is the projection operator and $A_{Q^1}:Q\cap C^1\rightarrow \ker\pi$ is the restriction of $\bar{A}.$ Using Taylor series expansion for $G(\Psi z+y, \nu)$ at $(z, y, \nu) = (0, 0, 0),$ above equation becomes

where

Now, we define the operator $N_j^1$ on $V_j^4(\mathbf{R}^2),$ the vector space of homogeneous polynomials in $z_1, \; z_2, \; \nu_1, \; \nu_2$ of degree $j$ having coefficients in $\mathbf{R}^2$, by

Using ^[18] and ^[20], we may decompose $V_j^4(\mathbf{R}^2)$ as follows:

Now, we consider a mapping from $V_j^4(\mathbf{R}^2)$ to $I_m(N_j^1)$ $(j = 2, 3)$ denoted by $P_{I, j}^1$ which satisfies

It is clearly visible that

Using ^[19,21], the normal form of Eq (2.2) on the center manifold associated to the space P can be taken as

where $f_2^1(z, 0, \nu) = Pr_{I_m(M_2^1)^c}g_2^1(z, 0, \nu)$

$= \left(\begin{array}{c} 0 \\ 2\phi_{22}\nu_1z_1+2[(\phi_{12}-\eta\phi_{22})\nu_1-\phi_{22}\nu_2]z_2+g{\prime\prime}(0)\phi_{22}z_1^2+2g{\prime\prime}(0)(\phi_{12}-\eta\phi_{22})z_1z_2 \\ \end{array} \right).$

The Eq (2.2) may be converted into the following normal form on the center manifold:

where $\lambda_1 = \phi_{22}\nu_1, \; \lambda_2 = (\phi_{12}-\eta\phi_{22})\nu_1-\phi_{22}\nu_2, \; c_{2} = \frac{g^{\prime\prime}(0)\phi_{22}}{2}, \; d_{2} = g^{\prime\prime}(0)(\phi_{12}-\eta\phi_{22}).$

We can have a small $a$ such that $\phi_{12}-\eta\phi_{22} < 0$ as $0 < a < M$. Without loss of generality, we are assuming that $\phi_{12}-\eta\phi_{22} < 0.$ As per our assumption $2+\epsilon(\eta-2\tau) > 0$ it is clear that $\phi_{22} > 0$. Hence, the sign of $g^{\prime\prime}(0)$ determines the sign of $c_2$ and $d_2$.

Case I: If $g^{''}(0) > 0,$ then $c_2 > 0$ and $d_2 < 0.$ After re-scaling the time and change the coordinates in the following manner:

The Eq (3.18) transforms to the following form on the center manifold

where $\alpha_1 = \frac{4(\phi_{12}-\eta\phi_{22})^2}{\phi_{22}}\nu_1$ and $\alpha_2 = -\frac{2(\phi_{12}-\eta\phi_{22})^2}{\phi_{22}}\nu_1+2(\phi_{12}-\eta\phi_{22})\nu_2$ (dropping bars).

The following summarizes the bifurcation curves related to the perturbation parameters $\nu_1, \; \nu_2$ ^[3,6,8,20]:

(i) If $\nu_1 = 0,$ then transcritical bifurcation occurs.

(ii) If $\nu_1 < 0$ and $3[2+\epsilon(\eta-2\tau)]\nu_2+[6\eta+\epsilon(3\tau^2+2\eta^2-6\tau\eta)]\nu_1 = 0$, then Hopf bifurcation from the zero equilibrium point ($H_0$) occurs.

(iii) If $\nu_1 > 0$ and $3[2+\epsilon(\eta-2\tau)]\nu_2-[6\eta+\epsilon(3\tau^2+2\eta^2-6\tau\eta)]\nu_1 = 0$, then Hopf bifurcation from the non-zero equilibrium point ($H_1$) occurs.

(iv) If $\nu_1 > 0$ and $21[2+\epsilon(\eta-2\tau)]\nu_2-5[6\eta+\epsilon(3\tau^2+2\eta^2-6\tau\eta)]\nu_1 = 0$, then Homoclinic bifurcation from the zero equilibrium point ($H_C^0$) occurs.

(v) If $\nu_1 < 0$ and $3[2+\epsilon(\eta-2\tau)]\nu_2+[6\eta+\epsilon(3\tau^2+2\eta^2-6\tau\eta)]\nu_1 = 0$, then Homoclinic bifurcation from the non-zero equilibrium point ($H_C^1$) occurs.

Case II: If $g^{\prime\prime}(0) < 0$ then $c_2 < 0$ and $d_2 > 0.$ After re-scaling the time and change the coordinates in the following manner:

The Eq (3.18) transforms to the following form on the center manifold (dropping the bars)

where $\alpha_1, \; \alpha_2$; and the corresponding bifurcation curves are the same as of Case I.

Case III: If $g^{\prime\prime}(0) = 0,$ then $c_2 = d_2 = 0$. We must compute higher order normal forms in order to talk about the dynamics near the B-T singularity. From ^[19] and ^[20], we have

where $\bar{g}_3^1(z, 0, \nu) = g_3^1(z, 0, \nu)+$$\frac{3}{2}[(D_zg_2^1)(z, 0, \nu)U_2^1(z, \nu)-$$(D_zU_2^1)(z, \nu)f_2^1(z, 0, \nu)+$$(D_yg_2^1)(z, 0, \nu)$$U_2^2(z, \nu)].$ After computation, we have

where $U_2^1$ and $U_2^2$ are defined in ^[20].

To compute $f_3^1(z, 0, \nu)$, we have used similar method as in our previous work ^[15].

Using the same process as in ^[20], we have the followings:

where $l_2^i(\theta)(z_1, z_2, \nu_1, \nu_2) = l_{22000}^i(\theta)z_1^2+l_{20200}^i(\theta)z_2^2+$$l_{20020}^i(\theta)\nu_1^2+l_{20002}^i(\theta)$$\nu_2^2 +l_{21100}^i(\theta)z_1z_2+l_{21010}^i(\theta)$$\nu_1z_1$ $+l_{21001}^i(\theta)\nu_2z_1+$$l_{20110}^i(\theta)\nu_1z_2+l_{20101}^i(\theta)$$\nu_2z_2+l_{20011}^i(\theta)\nu_1\nu_2$ for $i = 1, 2.$

From above equations, we have

where $\beta_1 = \phi_{22}\nu_1+\frac{1}{2}\phi_{22}l_{21010}^1(-\eta)\nu_1^2-\frac{1}{2}\phi_{22}l_{21010}^2(-\tau)\nu_1\nu_2$, $\beta_2 = (\phi_{12}-\eta\phi_{22})\nu_1-\phi_{22}\nu_2+$$\frac{1}{2}[\phi_{22}l_{20110}^1(-\eta)+\phi_{12}l_{21010}^1(-\eta)]\nu_1^2+$$\frac{1}{2}[\phi_{22}l_{20101}^1(-\eta)-\phi_{22}l_{20110}^2(-\tau)- \phi_{12}l_{21010}^2(-\tau)]$$\nu_1\nu_2-3\phi_{22}l_{20101}^2(-\tau)\nu_2^2,$ $c_3 = \frac{\phi_{22}}{6}[-6b+g^{\prime\prime\prime}(0)],$ $d_3 = \frac{1}{2}[\phi_{12}(-6b+g^{\prime\prime\prime}(0))-\eta\phi_{22}g^{\prime\prime\prime}(0)].$

Utilizing the following time re-scaling and co-ordinate shifting :

the Eq (3.24) converted to the following

where $\gamma_1 = (\frac{d_3}{c_3})^2\beta_1, \; \gamma_2 = -\frac{ d_3}{|c_3|}\beta_2, \; s = sgn(c_3).$

The sign of $s$ determines the bifurcations of the Eq (3.26) ^[22,23].

For s = 1, the bifurcations for the small perturbed parameters $(\nu_1, \, \nu_2)$ are summarised as follows ^[3,8,20]:

(i) On the curve $S = \{(\nu_1, \nu_2):\; \nu_1 = 0, \; \nu_2\in\mathbf{R}\},$ the Eq (3.26) possesses pitchfork bifurcation.

(ii) On the curve $H = \{(\nu_1, \nu_2):\; \nu_2 = -\frac{6\eta+\epsilon(3\tau^2+2\eta^2-6\tau\eta)}{3[2+\epsilon(\eta-2\tau)]}\nu_1+O(\nu_1^2), \; \nu_1 < 0\},$ the Eq (3.26) possesses Hopf bifurcation at the zero equilibrium point.

(iii) On the curve $L = \{(\nu_1, \nu_2):\; \nu_2 = -\frac{2[6\eta+\epsilon(3\tau^2+2\eta^2-6\tau\eta)]}{15[2+\epsilon(\eta-2\tau)]}\nu_1+O(\nu_1^2), \; \nu_1 < 0\},$ the Eq (3.26) experiences Heteroclinic bifurcation at the zero equilibrium point.

For $s = -1$, the bifurcation curves are concluded as follows ^[4,5,8]:

(i) On the curve $S = \{(\nu_1, \nu_2):\; \nu_1 = 0, \; \nu_2\in\mathbf{R}\},$ the Eq (3.26) experiences Pitchfork bifurcation.

(ii) On the curve $H_1 = \{(\nu_1, \nu_2):\; \nu_2 = -\frac{6\eta+\epsilon(3\tau^2+2\eta^2-6\tau\eta)}{3[2+\epsilon(\eta-2\tau)]}\nu_1+O(\nu_1^2), \; \nu_1 < 0\},$ the Eq (3.26) experiences Hopf bifurcation at the zero equilibrium point.

(iii) On the curve $H_2 = \{(\nu_1, \nu_2):\; \nu_2 = -\frac{4[6\eta+\epsilon(3\tau^2+2\eta^2-6\tau\eta)]}{3[2+\epsilon(\eta-2\tau)]}\nu_1+O(\nu_1^2), \; \nu_1 > 0\},$ the Eq (3.26) experiences Hopf bifurcation at the non-zero equilibrium point.

(iv) On the curve $T = \{(\nu_1, \nu_2):\; \nu_2 = -\frac{17[6\eta+\epsilon(3\tau^2+2\eta^2-6\tau\eta)]}{15[2+\epsilon(\eta-2\tau)]}\nu_1+O(\nu_1^2), \; \nu_1 > 0\},$ the Eq (3.26) experiences Homoclinic bifurcation.

(v) On the curve $H_d = \{(\nu_1, \nu_2):\; \nu_2 = -\frac{3.256[6\eta+\epsilon(3\tau^2+2\eta^2-6\tau\eta)]}{3[2+\epsilon(\eta-2\tau)]}\nu_1+O(\nu_1^2), \; \nu_1 > 0\},$ the Eq (3.26) experiences fold bifurcation of the limit cycle.

4.   Numerical results

We have only so far focused on the Eq (2.2)'s dynamical behaviour in terms of theoretical observations. This section discusses some numerical instances to support the validity of our theoretical conclusions and uncover other standout characteristics of the specified system. To show the phase diagrams, Mathematica 7.0, Maple 18, and MATLAB-R2015a were utilized as the computer programmes. The system parameter values are taken from the study done by Siewe et al. ^[24] with slightest deviation. Let us consider the positional delayed feedback function $g(x(t-\eta)) = \frac{\sin{x(t-\eta)}}{5},$ we have $g(0) = 0$, $g^{\prime}(0) = 0.2$, $g^{\prime\prime}(0) = 0.0$ and $g^{\prime\prime\prime}(0) = -0.2 < 0$. For the purpose of the simulation, we have selected the following set of parameters: $a = d = 0.2, \, b = 0.05, \, \epsilon = -0.1, \, \tau = 0.5$, $\eta = 0.5$, for which the condition of double zero eigenvalues (i) $\epsilon = -a\tau_2$ and (ii) $[2+\epsilon(\eta-2\tau)] = 2.05 > 0$ are fulfilled. Moreover, it is calculated that the value of $-\frac{6\eta+\epsilon(3\tau^2+2\eta^2-6\tau\eta)}{3[2+\epsilon(\eta-2\tau)]} = -0.4918699187$.

Following numerical studies display an agreement with our analytical findings: The phase portraits and the solution trajectory of Eq (2.2) with the critical B-T bifurcation parameter value $(d, \epsilon) = (0.2, -0.1)$ is demonstrated by Figure 1. Figure 1(a) corresponds to unperturbed Bogdanov-Takens bifurcation parameter $(d, \epsilon) = (0.2, -0.1)$ and Figure 1(b) is the phase diagram with a small perturbation $(\nu_1, \nu_2) = (-0.002, 0.00098)$. If we choose the initial point at $(0.06364, 0.06861)$, the system produce slowly converging diagram around the origin (see Figure 2).

We now take into account a minor alteration to the bifurcation parameters by letting $(d, \epsilon) = (0.2+\nu_1, -0.1+\nu_2)$. If we consider the perturbation parameters $(\nu_1, \nu_2) = (-0.01, 0.0)$, there exists a stable limit cycle around the trivial equilibrium of Eq (2.2) (see intersection of blue and red solution curves of Figure 3). From the bifurcation diagram (cf. Figure 4), one can be concluded that the the coefficient of damping $\epsilon$ has a crucial impact on regulating the oscillatory nature of Eq (2.2).

Next, we consider the parameter set $a = 0.02, \, b = 0.04, \, \epsilon = 0.002, \, \tau = 0.5, \, \eta = 0.0005$ and the feedback function $g(x) = \sin{x}$. A small perturbation $(\nu_1, \nu_2) = (-0.1, 0.1)$ to $(d, \epsilon)$ produce two homoclinic loop around the saddle point $(0, 0)$. In Figure 5, there exists a left homoclinic loop, where the trajectory started with the initial point $(-0.1, 0.1)$. The right homoclinic loop is also formed if the trajectory started with the initial point $(0.1, 0.1)$ (Figure is not reported here).

In Figure 6(a), (b) the largest Lyapunov exponents are depicted for the delayed feedback functions $\sin(x(t-\eta))$ and $\frac{e^{x(t-\eta)}-1}{e^{x(t-\eta)}+2}$ respectively for $\eta$ and other parameters remain the same as in Figure 5.

Figure 7 describes the Eq (2.2)'s chaotic dynamics. Figure 8 illustrates the existence of two small limit cycles covered by a larger limit cycle.

5.   Conclusions

In this present study, we have investigated the normal form for Bogdanov-Takens bifurcation experienced by the proposed retarded van der Pol-Duffing equation with multiple delay having possible degeneration. For the present system, we are aimed to study the dynamics with parametric positional delayed feedback in relation with quasi-periodic solution near and quasi-periodic attractor leading to chaotic solution. Analytically as well as numerically, we have found that the delayed Eq (2.2) experiences a B-T bifurcation if the conditions $\epsilon = -a\eta$ and $[2+\epsilon(\eta-2\tau)] > 0$ are satisfied simultaneosly. Utilizing the center manifold and normal form theories, we derived the canonical forms of Bogdanov-Takens singularity. Moreover, the phase and bifurcation diagrams for different positional delayed feedback functions are shown. It is shown that that the the coefficient of damping plays an important function to regulate dynamics of the Eq (2.2) (cf. Figure 4). One of the most crucial mechanisms for examining the chaotic behaviour of dynamical systems is the Lyapunov exponent. The positivity of the computed values of Lyapunov exponents provide a definite conclusion of chaotic dynamics. However, Lyapunov exponents are also notoriously difficult to estimate reliable from experimental data. One should be cognizant of additional evidence proving Sensitive Dependence on Initial Conditions (SDIC) before asserting that any system is chaotic. Growing awareness of this challenge will undoubtedly continue to encourage the creation of newer and more reliable algorithms for finding the Largest Lyapunov Exponents. From the Figure 6(a), (b), it is seen that the Largest Lyapunov exponents are at the positive level for the positional feedback function $\sin(x(t-\eta))$ and $\frac{e^{x(t-\eta)}-1}{e^{x(t-\eta)}+2}$ respectively, when $\eta$ is increased to $20.5$ (other parameters remain the same as in Figure 5). For the both cases, similar chaotic dynamics have been observed and shown in Figure 7. Figure 8 shows three close limit cycles encircling the trivial equilibrium point and other two symmetrically situated axial equilibria corresponding to very neighboring starting position, which exists only for the end of the total time span of evolution of the trajectories. We employ the MatLab version of the approach developed by Wolf et al. ^[25] to estimate the largest Lyapunov exponent for the current study. Multiple delays, as is well known, make it difficult and tiresome to analyze the distribution of the eigenvalues and cause the system to exhibit considerably deeper dynamical characteristics. There are hardly any unfolding results about the double Hopf bifurcation, quasi-periodic attractors and coexisting attractors enveloping. The engineers can use these results as a guide for selecting the delay settings to obtain the desired dynamical effects. Therefore, there are less number of results about the quasi-periodic solution leading to chaos due to multiple delays for dynamical systems on oscillators. We'll work about these realities in the next study.

Acknowledgements

The authors gratefully acknowledge Qassim University, represented by the Deanship of Scientific Research, on the financial support for this research under the number 10367-cos-2020-1-3-I during the academic year 1442 AH / 2020 AD.

Conflict of interest

The authors declare no conflicts of interest.