# Bifurcation structure of two coupled periodically driven double-well Duffing oscillators

Anatole Kenfack

► To cite this version:

Anatole Kenfack. Bifurcation structure of two coupled periodically driven double-well Duffing oscillators. Chaos, Solitons & Fractals, 2003, 15 (2), pp.205-218. 10.1016/S0960-0779(01)00250-8. hal-01356790

HAL Id: hal-01356790

https://hal.science/hal-01356790

Submitted on 26 Aug 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Bifurcation structure of two coupled periodically driven double-well Duffing oscillators

Anatole Kenfack

Max-Planck-Institute for the Physics of Complex Systems, Nöthnitzer Strasse 38, D-01187 Dresden, Germany
Department of Physics, University of Dschang, P.O. Box 67, Dschang, Cameroon

The bifurcation structure of coupled periodically driven double-well Duffing oscillators is investigated as a function of the strength of the driving force $f$ and its frequency $\Omega$. We first examine the stability of the steady-state in linear response, and classify the different types of bifurcation likely to occur in this model. We then explore the complex behavior associated with these bifurcations numerically. Our results show many striking departures from the behavior of coupled driven Duffing oscillators with single-well potentials, as characterized by Kozłowski et al. [Phys. Rev. E 51 (1995) 1861]. In addition to the well-known routes to chaos already encountered in a one-dimensional Duffing oscil-lator, our model exhibits imbricated period-doubling of both types, sym m etry-breaking,sudden chaos and a great abundance of Hopf bifurcations, many of which occur more than once for a given driving frequency. We explore the chaotic behavior of our model using two indicators, namely Lyapunov exponents and the power spectrum. Poincaré cross-sections and phase portraits are also plotted to show the m anifestation of coexisting periodic and chaotic at-tractors including the destruction of $T^2$ tori doubling.

## 1. Introduction

In recent years, a large number of theoretical calculations, numerical simulations and experiments have been carried out on systems of coupled anharmonic oscillators which provide fundamental models of the dynamical problems in several disciplines. Such models describe a wide range of processes, helping in general to understand the routes to chaos that take place in biological, chemical, physical, electromechanical and electronical systems. Here we refer only to recent papers close to the subject of our own [1–18]. Coupling may arise between oscillators of the same types or oscillators of different types. Amongst these systems, the most intensively investigated examples are the Duffing oscillators [1–8] (strictly dissipative systems which converge to quiescence when not driven), the Van der Pol oscillators [9–11] (self-excited systems whose final behavior when not driven is a limit cycle) and the coupled Van der Pol–Duffing oscillators [12,13].

In this paper, we present the results of an investigation of two identical coupled double-well Duffing oscillators subjected to a periodically driven force. Our study was stimulated by the earlier work of Kozłowski et al. [1] who have studied the model of the same type, but considered only single-well potentials. Using bifurcation diagrams and phase diagrams these authors have demonstrated that the global pattern of bifurcation curves in parameter space consists of repeated subpatterns similar to the superstructure found in a one-dimensional driven Duffing oscillator. They have also
proven the existence of a Hopf bifurcation which does not appear in a model of one-dimensional Duffing oscillator. Other interesting results have been reported in the context of coupled Duffing oscillators. In [2] Kunick et al. showed that two coupled oscillators can behave regularly with nonzero coupling and chaotically with zero coupling. In a model of two coupled Duffing oscillators driven with incommensurate frequencies and coupled additively, Stagliano et al. [3] observed the period-doubling of an attracting T² torus and its destruction in parameter space. Migration control in two Duffing oscillators by open-plus-closed-loop control method and adaptive control algorithm were studied by Paul Raj and Rajasekar [4]. These authors, in another system of two coupled Duffing oscillators, have also investigated a basin of attraction with several coexisting chaotic attractors and synchronization of chaos [5]. Yagasaki [6], after having described an Auto driver Hommap for the numerical analysis of homoclinic in maps of periodically forced Duffing systems, modified a version of the homoclinic Melnikov method for orbits homoclinic to two types of periodic orbits. These theories have been successfully applied to weakly coupled Duffing oscillators. Moreover Yin et al. [7] investigated the effect of phase difference on the driving forces of two coupled identical Duffing oscillators. They observed desynchronization and lag synchronization of chaotic attractors. Bifurcation behaviors, showing a kind of Hopf bifurcation, from synchronous chaos of a chain of Duffing oscillators, have been explored by Ma Wen-Qi et al. [8] using the generalized winding number in tangent space.

The equations of motion for the two identical coupled double-well Duffing oscillators that we are interested in are the following:

$$ \frac{d^2 x}{dt^2} = -\alpha \frac{dx}{dt} + x - x^3 + k(y - x) + f \cos(\Omega t), \quad (1) $$

$$ \frac{d^2 y}{dt^2} = -\alpha \frac{dy}{dt} + y - y^3 - k(y - x), $$

where $\alpha$ is the damping parameter, $k$ the coupling parameter, $f$ and $\Omega$ are the amplitude and the frequency of the driving force, respectively. In system (1), we represent the oscillators by the state variables $(x, dx/dt)$ as oscillator A or sub-

<figure><img src="image_2.png" /><figcaption>Fig. 1. Eigenvalues spectrum for the stability of the system with $\alpha \in [-0.2, 0.2]$ and $k \in [-5, 5]$. Note that when $\alpha < 0$, the eigenvalues are in the half-plane left and in the half-plane right otherwise.</figcaption></figure>
system A subjected to the periodic force $f(t) = f \cos(\Omega t)$, which is coupled to the other one with the state variables $(y, dy/dt)$ as oscillator B or subsystem B. Here $f$ and $\Omega$ are the control parameters of our system. As evident, for $k=0$, Eq. (1) describes two uncoupled Duffing oscillators with their motion governed by $\alpha$ and $f$. The coupling considered can be interpreted as a perturbation of each oscillator through a signal proportional to the difference of their amplitudes. A natural question to ask is how chaotic attractors arise as these parameters are varied. In the case of a one-dimensional double-well Duffing oscillator, quasiperiodic and period-doubling routes to chaos have been found [19–27].

<figure><img src="image_3.png" /><figcaption>Fig. 2. Bifurcation diagrams showing the first coordinate $x$ for the driving amplitude $f = 15$ (a) and the second coordinate $v$ for $f = 20$ (b) of the Poincaré section vs the driving frequency $\Omega$. Windows of periodic solutions and chaotic domains are visible. Near $\Omega = 0.385$ (a), six pairs of *sh* occur followed by the reversed *pd* till the chaotic domain. Imbricated *pd*, *sh*, *sn* are mentioned. Note that similar behavior as (b) is observed when $\Omega \le 0.4$ and for small values of $f$ ($f < 15$).</figcaption></figure>
The rest of the paper is structured as follows. In Section 2, we establish a linear stability analysis of the system projected on the Poincaré map. Section 3 is devoted to the survey of bifurcation diagrams and the chaotic behavior is characterized in Section 4. We end with the conclusions of the work in Section 5.

## 2. Stability analysis

For stability analysis and even for purposes of numerical simulations, it is convenient to transform the second-order differential Eq. (1) into an autonomous system of first-order differential equations of the following form:

<figure><img src="image_4.png" /><figcaption>Fig. 3. Bifurcation diagrams for the driving force $f = 4.5$ (f), $f = 15$ (a,b,e), $f = 20$ (d), $f = 25$ (c) with $\Omega > 0.4$. Throughout imbricated *sb*, *sn*, *pd* of both types and chaotic domains are observed. Besides, resonances (R) appear at $\Omega \approx 0.475$ (a,b), $\Omega \approx 0.445$ (d), $\Omega \approx 0.555$ (e) and sudden chaos at $\Omega \approx 0.408$ (c).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 3 (continued)</figcaption></figure>

$$
\begin{align*}
\frac{dx}{dt} &= v, \\
\frac{dv}{dt} &= -\alpha v + x - x^3 + k(y - x) + f \cos(\theta), \\
\frac{dy}{dt} &= w, \\
\frac{dw}{dt} &= -\alpha w + y - y^3 - k(y - x), \\
\frac{d\theta}{dt} &= \frac{\Omega}{2\pi}
\end{align*}
\tag{2}
$$

or equivalently

$$
\frac{dV}{dt} = F(V, \mu),
$$

where $T = 2\pi/\Omega$ is the period of the driving force, $\theta$ the cyclic variable, $V(x, v, y, w, \theta)$ an autonomous vector field and $\mu(\alpha, k, f, \Omega)$ is an element of the parameter space. This system generates a flow $\phi = \{\phi^T\}$ on the phase space $\mathbf{R}^4 \times \mathbf{S}^1$ and there exists a global map:

$$
\begin{align*}
P: \quad \Sigma_c &\rightarrow \Sigma_c, \\
V_p(x, v, y, w) &\rightarrow P(V_p) = \{\phi^T\}_{\Sigma_c(x,v,y,w,\theta_0)}
\end{align*}
$$

with $\theta_0$ being a constant determining the location of the Poincaré cross-section and $(x, v, y, w)$ the coordinates of the attractors in the Poincaré cross-section $\Sigma_c$ defined by:

$$
\Sigma_c = \{(x, v, y, w, \theta) \in \mathbf{R}^4 \times \mathbf{S}^1 | \theta = \theta_0\}.
$$

System (2) is symmetric since the transformation

$$
S : (x, v, y, w, \theta) \rightarrow (-x, -v, -y, -w, \theta + \pi)
$$
leaves (2) invariant. This system can support symmetric orbits and also asymmetric ones which are not invariant with the transformation S. We use a perturbation analysis to proceed with the analysis of solution stability in the Poincaré map $\Sigma_c$. This method consists of differentiation and translation. Thus if we add a small perturbation $\delta V_p(\delta x, \delta v, \delta y, \delta w)$ to the steady-state $V_{p0}(x_0, v_0, y_0, w_0)$ we obtain equations for perturbed motion

$$x = x_0 + \delta x, \quad y = y_0 + \delta y, \quad v = v_0 + \delta v, \quad w = w_0 + \delta w. \tag{3}$$

Substituting Eq. (3) into Eq. (2) and developing into powers of the above small perturbations yield the following variational differential equation:

$$\frac{d\delta V_p}{dt} = DF(V_{p0})\delta V_p \tag{4}$$

with

$$DF(V_{p0}) = \begin{pmatrix} 0 & 1 & 0 & 0 \\ 1 - k - 3x_0^2 & -\alpha & k & 0 \\ 0 & 0 & 0 & 1 \\ k & 0 & 1 - k - 3y_0^2 & -\alpha \end{pmatrix},$$

where $DF(V_{p0})$ is the Jacobian $4 \times 4$ matrix describing the vector field along the solution $\delta V_p(t)$ with $V_{p0}$ being an equilibrium point or the steady-state. Next for the very small amplitudes, we neglect in the matrix $DF(V_{p0})$ the powers higher than one, that is, $x_0^2$ and $y_0^2$. The solution after one period $T$ of the oscillations in the linearized Poincaré map is simply expressed as:

$$\delta V_p(T) = \delta V_{p0} \exp(DF(V_{p0})T), \quad \delta V_p(0) = \delta V_{p0}, \tag{5}$$

where $DF(V_{p0})$ is the monodromy matrix of a periodic orbit connecting arbitrary infinitesimal variations in the initial conditions $\delta V_{p0}$ with corresponding changes $\delta V_p(T)$, after one period $T$. The stability of the periodic motion is therefore determined according to the real parts of the roots of the characteristic equation $\det(DF(V_{p0}) - I\lambda) = 0$ written in the following form:

$$\lambda^4 + A_3\lambda^3 + A_2\lambda^2 + A_1\lambda + A_0 = 0, \tag{6}$$

where $(A_i)$, $i = 0, 1, 2, 3$ are the coefficients depending on the two parameters $k$, $\alpha$ and $\lambda = (\lambda_i)$ the eigenvalues of $DF(V_{p0})$. The roots of Eq. (6) are obtained using the Bairstow–Newton–Raphson algorithm [28] and for different values of $k \in [-5, 5]$ and $\alpha \in [-0.2, 0.2]$. The results as the spectrum of the eigenvalues are depicted in a complex plane $C$ (see Fig. 1). This figure enables us to get an idea of the different types of bifurcation likely to appear in the system. Since $DF(V_{p0})$ is a real matrix, complex eigenvalues occur in a complex conjugate pairs responsible for the symmetry observed along the real axis. Thus if $\lambda_i$ is real, it is clear from Eq. (5) that the eigenvalue is the rate of contraction $(\lambda_i < 0)$ or expansion $(\lambda_i > 0)$ near the steady-state. Next if $\lambda_i$ is complex, the real part of $\lambda_i$ gives the rate of contraction $(\mathrm{Re}(\lambda_i) < 0)$ or expansion $(\mathrm{Re}(\lambda_i) > 0)$ of the spiral while the imaginary part $\mathrm{Im}(\lambda_i)$ is of the frequency rotation. This can be well understood with the following expression of the eigenvalues of the linearized Poincaré map:

$$\sigma_i = \exp(\mathrm{Re}(\lambda_i)) [\cos(\mathrm{Im}(\lambda_i)T) + j \sin(\mathrm{Im}(\lambda_i)T)]. \tag{7}$$

Globally, if $\mathrm{Re}(\lambda_i) < 0$ for all $\lambda_i$, then all sufficiently small perturbations tend toward zero as $t$ goes to infinity and the steady-state (nodes $(n)$, saddle-nodes $(sn)$, spiral $(sp)$) is stable. If $\mathrm{Re}(\lambda_i) > 0$ for all $\lambda_i$, then any small perturbation grows with time and the steady-state $(n, sn, sp)$ is unstable. (A stable or unstable equilibrium state with no complex eigenvalues is often called node.) On the other hand if there exist $i$ and $l$ such that $\mathrm{Re}(\lambda_i) < 0$ and $\mathrm{Re}(\lambda_l) > 0$, therefore the equilibrium state is nonstable. (A nonstable equilibrium state is often called saddle; an equilibrium point whose eigenvalues all have a nonzero real part are called hyperbolic [29].) With regard to the precedent with again a look on the spectrum of Fig. 1, it follows that our system of the given parameter $k$ and $\alpha$ can undergo many types of bifurcations namely: saddle-node ($\mathrm{Re}(\lambda_i) = 1$), period-doubling ($\lambda_i = -1$), Hopf bifurcations ($\lambda_i = \gamma \pm j\beta$ ($j^2 = -1$) with $\gamma < 0$ and the motion has the form of exponential increasing beats with period $2\pi/\beta$ [17]) and of course symmetry-breaking bifurcation which is often a prerequisite for the first period-doubling bifurcation [30]. Except Hopf bifurcation, such bifurcations have been successfully found in one-dimensional double-well Duffing oscillators subjected to a periodically driven force [19,21,26] using the eigenvalues called Floquet multipliers of the linearized Poincaré map. In a coupled identical single-well Duffing oscillators, Kozłowski et al. [1] showed by linear analysis a similar eigenvalues scenario and they conjectured that the resulting alternating bifurcation sequences have to be expected for all systems of coupled strictly dissipative oscillators that are driven periodically. The above analysis allows us to know how the steady-state
becomes locally unstable and to be aware of the type of bifurcation expected in the system. Here, no exploding amplitude is possible since the shape of the potential offers globally bounded solutions ($V(x, y) = -(1/2)x^2 - (1/2)y^2 + (1/4)x^4 + (1/4)y^4 \to +\infty$ as $|x|, |y| \to \infty$).

### 3. Local bifurcations

The aim of this section is to seek numerically the routes which lead to chaotic solutions of the system when the control parameter $\Omega$ evolves, for different values of the driving force $f$. Bifurcation diagrams are very good for numerical as well as for experimental studies when there is a tunable parameter. When a control parameter is varied and bifurcation takes place a qualitative change of the system happens. For the numerical computations of these diagrams the control parameter $\Omega$ is increased from an initial value $\Omega_i$ to a final value $\Omega_f$ and then decreased from $\Omega_f$ to $\Omega_i$ in a very small step ($\Delta\Omega = 10^{-5}$). The last computed cyclic point for a given value of $\Omega$ is always used as a new initial value for the next value of $\Omega$. Starting with initial conditions $(x, v) = (1, 0)$, $(y, w) = (1, 0)$ at $\Omega_i$, system (2) is integrated, using the standard Runge-Kutta Algorithm [28], for 100 periods of the driving force until the transient has died out; the trajectories expected are attractors and the local calculation error is sufficiently small. Then to find out whether the trajectory is periodic (quasiperiodic) or chaotic, the system is integrated for the next 180 periods in order to catch the maximum of coexisting attractors. This procedure allows us to make sure that coexisting attractors are realized due to different initial conditions which are handed throughout from one parameter to another. We chose $k = 5$ and $\alpha = 0.1$ for all examples in the subsequent section. The bifurcation diagrams obtained here show the projection of attractors in the Poincaré section onto $x$ or $v$ vs the control parameter $\Omega$ and for different values of $f$. In a one-dimensional double-well Duffing oscillator [19–21,23–27], it is shown that chaotic solutions arise through quasiperiodic, period-doubling cascade with a sequence of symmetry-breaking and saddle-nodes. Hopf bifurcations have not been found, but are abundant in our model.

We computed a great number of diagrams, and some of them showing the main features of the system are presented in Figs. 2–4. It is important to note that no hysteresis has been detected while checking all these diagrams. We can therefore avoid any confusion over the periods of attractors when reading this type of diagram; for instance, one single period-2 yields to two points for one frequency value whereas two coexisting attractors each of period-1 also yield to two points for one frequency. These diagrams show a great number of coexisting attractors (chaotic domains) intermingled with imbricated windows made up of periodic solutions of different periodicities, period-doubling of both types, sudden Chaos and Hopf bifurcation. The so-called chaotic domains (or chaotic sea) will be characterized in the following section.

At lowest frequencies ($0 < \Omega \le 0.4$) quasiperiodic as well as chaotic solutions occur together with a number of symmetry-breaking ($sb$), period-doubling ($pd$) and saddle-nodes ($sn$). In Fig. 2(a) where $f = 15$, it appears in the neighborhood of $\Omega = 0.385$, six pairs of $sb$ followed by three reversed $pd$ until the chaotic domain. When the driving force is increased to $f = 20$, overlapping $pd$ cascade and $sn$ appear as can be seen in Fig. 2(b). Similar features as in the latter case are also observed for small values of $f$ ($f < 15$). In other ranges of driving frequencies bigger than $\Omega = 0.4$ the aforementioned features still exist and additional ones such as resonances and Hopf bifurcations appear in the system; the results are displayed in Figs. 3 and 4. In Fig. 3 it can be seen throughout overlapping $pd$ cascade of both types, $sb$, resonances ($R$) and even sudden chaos ($SC$). When $f = 4.5$ in Fig. 3(f) at $\Omega \approx 0.785$, period-8 attractor undergoes reversed period-doubling cascade and period-2 is created as from $\Omega \approx 0.805$. Besides, for $f \ge 7$ appear resonances ($R$) in Figs. 3(a), (b), (d), and (e), Fig. 6(a) and sudden chaos in Fig. 3(c). However in Fig. 4 when $f$ becomes much larger and $\Omega \ge 0.5$, Hopf bifurcation ($H$) is abundant and often appears more than once at a given frequency. This notation – $nH$– means that there are $n$ ($n \le 3$), integer Hopf bifurcation at a given frequency. One can therefore observe –$H$– in (a), $-3H$– in (b), $-2H$– $-2H$– in (c,d), $-2H$– $-3H$– $-2H$– in (e,f) and $-3H$–$3H$– in (g,h). In Fig. 4(h), for instance, at $\Omega = 0.5925$, a period-6 attractor undergoes reversed period-doubling to become a period-3; next after three pairs of symmetry-breaking near $\Omega \approx 0.6$, appears again period-3 attractor as from $\Omega \approx 0.601$. And then at $\Omega = 0.603$ this attractor of period-3 undergoes three Hopf bifurcations. Windows of periodic solutions separated by quasiperiodic ones are clearly visible. The Poincaré cross-section of quasiperiodic attractors of this kind ($-3H$–) consists of three invariant tori (destroyed or not). Thus at any $\Omega$ where $-nH$– exists, the projection of the attractor corresponds to $n$ Tori on the Poincaré cross-section. The torus attractor arises from Hopf bifurcation which is a bifurcation from a fixed point to an invariant curve. At the transition to quasiperiodic motion, this fixed point (projection of the limit cycle into the cross-section) loses its stability and gives birth to an invariant circle, which is a cross-section of the torus in the flow. The three pairs of $sb$ appearing at the neighborhood of $\Omega \approx 0.6$ (Figs. 4(g) and (h)) seem to be apparently a window of period-6 attractor. In fact it is just a period-3 attractor as can be illustrated further in Fig. 8(b). This orbit, born asymmetrically, coexists with its inversion symmetric orbit of the same period-3 and are caught simultaneously.
<figure><img src="image_6.png" /><figcaption>Fig. 4. Bifurcation diagrams for *f* = 30 (a,g,h), *f* = 25 (b), *f* = 15 (c,d), *f* = 20 (e,f) showing many Hopf bifurcations (*H*). One can observe *−H*− (a), *−3H*− (b), *−2H*−...−2*H* (c,d), *−2H*−...−3*H*−...−2*H*− (e,f), *−3H*−3*H*− (g,h) (see also text).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 4 (continued)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5. Bifurcation diagram for $f = 5$ showing attractors of small number of periods. Note that the system becomes regular and this behavior is observed for $\Omega \ge 3$ and for any value of $f$.</figcaption></figure>

Next we pursue our investigation for large values of $\Omega$ ($\Omega \ge 3$) and for different values of $f$. It turns out that the features observed so far have disappeared and the system becomes regular providing attractors with a small number of periods as can be seen in Fig. 5. The phase diagrams which often show bifurcations, curves, surfaces and which is out of the scope of this paper could also help to match the features obtained in this analysis.

## 4. Chaotic behavior

To better support the results obtained above, we intend in this section to characterize chaotic behavior. For this purpose some indicators are used namely the Lyapunov exponents and the power spectrum density. Also, Poincaré
cross-sections and phase portraits are shown to illustrate some attractors present in the system. The Lyapunov exponent is one of the most important tools for understanding chaotic behavior, obtained by examining “a very sensitive dependence of initial conditions”. In particular it is generally well established that a rigorous measure of chaos may be given in terms of Lyapunov spectrum of the dynamical system. A positive Lyapunov exponent is characteristic of chaos while zero and negative values of the exponent signify a marginally stable or quasiperiodic orbit and periodic orbit, respectively. Solving numerically the variational equation (4) together with system (2) we calculate the maximum Lyapunov exponent $\lambda_{max}$ in the Poincaré cross-section with

<figure><img src="image_9.png" /><figcaption>Fig. 6. Bifurcation diagram for $f = 15$ (a) and the corresponding Lyapunov spectrum (b) in the Poincaré map. Windows of chaotic domains are clearly justified with positive values of $\lambda_{max}$.</figcaption></figure>
$$ \lambda_{\max} = \lim_{\tau \to +\infty} \frac{1}{\tau} \ln(\|L(\tau)\|), $$

where $\|L(\tau)\| = (\delta x^2 + \delta v^2 + \delta y^2 + \delta w^2)^{1/2}$. The spectrum is displayed in Fig. 6 together with a bifurcation diagram. Windows of chaotic domains and periodic solutions are clearly justified. Power spectrum density using fast fourier transformation (FFT) has also been computed with some values from Fig. 6 for a periodic solution ($\Omega = 0.495, f = 7$) and a chaotic solution ($\Omega = 0.47, f = 7$) (see Fig. 7). Figs. 8 shows the phase portraits and Poincaré cross-section. In Figs. 8(a) and (b) plotted are the phase portraits, related to Figs. 4(g) and (h), showing a period-3 attractor for $f = 30$, $\Omega = 0.6025$ (a) and a coexistence of two asymmetric attractors, each of period-3 for $f = 30$, $\Omega = 0.6$ (b). The dots (●)

<figure><img src="image_10.png" /><figcaption>Fig. 7. Power spectra densities (PSD). The values of $f = 7$ and $\Omega$ are from the bifurcation diagrams shown in Fig. 6 with $\Omega = 0.495$ (a) periodic motion and $\Omega = 0.47$ (b) chaotic motion.</figcaption></figure>
represent the corresponding points in the Poincaré cross-sections. In Figs. 8(c) and (d) two chaotic attractors are plotted in the Poincaré cross-section for $\Omega = 0.5$, $f = 2$ and for $\Omega = 0.375$, $f = 15$, respectively. Furthermore, we observe in Figs. 8(e) and (f) three $T^2$ tori doubling at $\Omega = 0.605$, $f = 30$ and at $\Omega = 1.14$, $f = 15$, two $T^2$ tori doubling (destroyed) corresponding to $-3H$- and $-2H$-, respectively (values from Figs. 4(c) and (h)). In the torus route to chaos, the original torus appears to split into two circles at the torus doubling bifurcation point. This route is reminiscent of the period-doubling route to chaos, although there are a finite number of torus doublings before the onset of chaos motion. Also, the torus doubling route to chaos is a higher-dimensional phenomenon, requiring at least a four-dimensional flow or a three-dimensional map. It is not observed in one-dimensional maps, unlike the period-doubling route to chaos [30].

<figure><img src="image_11.png" /><figcaption>Fig. 8. Phase portraits and projection of the attractors in the Poincaré section onto the two first coordinates $x, v$ of the Poincaré cross-section. At $\Omega = 0.6025$ (a), one period-3 orbit and at $\Omega = 0.6$ (b) coexistence of two asymmetric orbits of the same period-3. The values of $\Omega$ are from Figs. 4(g) and (h) and $f = 30$. The dots (●) represent the corresponding points in the Poincaré cross-section. Two chaotic attractors shown in (c,d) for $\Omega = 0.5$, $f = 2$ and $\Omega = 0.375$, $f = 15$, respectively. Three $T^2$ tori doubling for $\Omega = 0.605$, $f = 30$ (e) and two $T^2$ tori doubling (destroyed) for $\Omega = 1.14$, $f = 15$ (f). The values of $\Omega$ are from Figs. 4(c) and (h).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 8 (continued)</figcaption></figure>

**5. Conclusion**

To sum up we have investigated coupled double-well Duffing oscillators subjected to a periodic driving force. Our main conclusion is that the shape of the potential in the coupled Duffing oscillators system has a profound influence on the routes to chaos taken by the system. By a linear analysis we have examined the stability of the steady-state of the system leading to different types of bifurcation likely to appear for $k \in [-5, 5]$ and $\alpha \in [-0.2, 0.2]$. With the help of bifurcation diagrams, it has been shown that the bifurcation structure, which is complicated, depends strongly on the values of the control parameter $\Omega$. For a simple electronic realization of our model, it is obvious that the oscillators may be very difficult to control for small values of $\Omega$ due to the large windows for chaos. This system becomes regular for large values of $\Omega$ ($\Omega \ge 3$). Apart from the routes already encountered in a one-dimensional double-well Duffing oscillator such as quasiperiodic and period-doubling cascade, sudden chaos and mostly Hopf bifurcations have been newly found in this work. Besides, resonances ($R$), imbricated $sh$, $sn$, $pd$ intermingled by chaotic domains are observed, rending the structure highly chaotic. Two and three $T^2$ tori doubling have also been observed. To characterize chaotic behavior of this system, Lyapunov exponents and the power spectrum density using FFT have been employed. Poincaré cross-sections and phase portraits have also been used for a better visualization of the periodic and chaotic attractors. An analytical approach using approximative technics could be an interesting topic for future work, and would probably shed further light on the rich behavior of this model.

**Acknowledgements**

It is a pleasure to thank Prof. Dr. Jan M. Rost for the illuminating discussions and Dr. Nic Shannon for a careful reading of the manuscript. Financial support by Alexander von Humboldt (AvH) Foundation/Bonn-Germany, under the grant of Research fellowship No.IV.4-KAM 1068533 STP, is gratefully acknowledged.

**References**

[1] Kozłowski J, Parlitz U, Lauterborn W. Phys Rev E 1995;51:1861.

[2] Kunick A, Steeb WH. J Phys Soc Jpn 1985;54:1220.

[3] Stagliano Jr. JJ, Wersinger JM, Slaminka EE. Phys D 1996;92:164.
[4] Paul Raj S, Rajasekar S. Phys Rev E 1997;55:6237.

[5] Paul Raj S, Rajasekar S, Murali K. Phys Lett A 1999;264:283.

[6] Yagasaki K. Int J Bifur Chaos 1998;8(7):1617;
Nonlinear Dyn 1999;20:319.

[7] Yin H-W, Dai J-H. Phys Rev E 1998;58:5683.

[8] Wen-Qi Ma, Jun-Zhong Y, Wen-Ji L, Gang B, Gang H. Acta Phys Sinica 1999;48(5):787.

[9] Poliashenko M, Mckay SR, Smith CW. Phys Rev A 1991;43:5638;
Phys Rev A 1991;44:3452;
Phys Rev A 1992;46:5271.

[10] Kapitaniak T, Steeb WH. Phys Lett 1991;152:33.

[11] Pastor I, Pérez-García VM, Encinas-Sanz F, Guerra JM. Phys Rev E 1993;48:171;
Phys Rev E 1995;52:1480.

[12] Woafo P, Chedjou JC, Fotsin HB. Phys Rev E 1996;54:5929.

[13] Han Y-J. J Korean Phys Soc 2000;37(1):3.

[14] Nath A, Ray DS. Phys Rev A 1986;34:4472.

[15] Umberger DK, Gregori C, Ott E, Afeyan B. Phys Rev A 1989;39:4835.

[16] Dressler U, Lauterborn W. Phys Rev A 1990;41:6702;
Dressler U. Phys Rev A 1998;38:2103.

[17] Püst L, Szöllös O. Int J Bifur Chaos 1999;9:479.

[18] Matthews PC, Mirollo RE, Strogatz SH. Phys D 1991;52:293.

[19] Englisch V, Lauterborn W. Phys Rev A 1991;44:916.

[20] Pecora L, Carroll TL. Phys Rev Lett 1991;67:945.

[21] Wang CS, Kao YH, Huang JC, Gou YS. Phys Rev A 1992;45:3471.

[22] Zeni AR, Gallas JAC. Phys D 1995;89:71.

[23] Osipov G, Glatz L, Troger H. Chaos, Solitons and Fractals 1998;9:307.

[24] Litvak-Hinenzon A, Rom-Kedar V. Phys Rev E 1997;55:4964.

[25] Chacón R, García-Hoz AM. Phys Rev E 1999;59:6558.

[26] Kim S-Y. Int J Mod Phys B 2000;14(17):1801;
Kim S-Y, Kim Y. Phys Rev E 2000;61:6517.

[27] Szemplińska-Stupnicka W, Zubrzycki A, Tyrkiel E. Int J Bifur Chaos 2000;10(6):1367.

[28] Press WH, Teukolsky SA, Vetterling WT, Flannery BP. In: Numerical recipes in Fortran 77. Cambridge: Cambridge University
Press; 1992. p. 370.

[29] Parker TS, Chau LO. In: Practical numerical algorithms for chaotic systems. New York: Springer; 1992. p. 57–82.

[30] Tufillaro NB, Abbot T, Reilly J. An experimental approach to nonlinear chaos. Reading, MA: Addison-Wesley; 1992.