We explore the dynamics of a periodically driven Duffing resonator coupled elastically to a van der Pol oscillator in the case of 1 : 1 internal resonance in the cases of weak and strong coupling. Whilst strong coupling leads to dominating synchronization, the weak coupling case leads to a multitude of complex behaviours. A two-time scales method is used to obtain the frequency-amplitude modulation. The internal resonance leads to an antiresonance response of the Duffing resonator and a stagnant response (a small shoulder in the curve) of the van der Pol oscillator. The stability of the dynamic motions is also analyzed. The coupled system shows a hysteretic response pattern and symmetry-breaking facets. Chaotic behaviour of the coupled system is also observed and the dependence of the system dynamics on the parameters are also studied using bifurcation analysis.

1. Introduction

In the macro world, one can easily find models of single, and coupled nonlinear oscillators in many different scientific disciplines, for example, biology, electronics and physics. However, in the small world, recent development in M/NEMS (Micro/Nano Electromechanical Systems) technology now can enable physicists to design and manufacture large arrays of coupled M/NEMS nonlinear oscillators to explore the mathematics and physics of emergent phenomena that appear in the self organization of the interacting individual elements. For instance, intrinsic localized modes and moving discrete breathers have been observed in coupled microcantilever resonators [1, 2]. Therefore, the fundamental study of coupled nonlinear oscillators is significant in understanding the emergent behaviour of complex dynamical systems and developing novel M/NEMS devices [35]. This work is also motivated by the fact that analysis of simpler cases as the building blocks can help us to gain insight into larger complicated systems [6, 7]. Among the building blocks studied in the literature, the essential elements are either self-sustained oscillators (van der Pol oscillators) or dissipative oscillators (Duffing type resonators) and the most intensively studied cases are the coupled van der Pol and the coupled Duffing oscillators [810]. To the best of our knowledge, less has been done in a dynamic system consisting of a Duffing type resonator coupled to a self-excited oscillator. Our aim in this paper is to consider the dynamics of such a system externally driven by a periodic force. The motion of the coupled dynamical system in dimensionless form is described by the following set of equations: ̈𝑥1+𝑐1̇𝑥1+𝑥1+𝑘3𝑥31=𝑘𝑐𝑥2𝑥1+𝑓𝑑𝜔cos𝑑𝜔1𝜏,̈𝑥2+𝑐2𝑥221̇𝑥2+𝜔2𝜔12𝑥2=𝑘𝑐𝑥1𝑥2,(1.1) where 𝑥1, 𝑐1, 𝜔1 and 𝑥2, 𝑐2, 𝜔2 are the displacement, damping coefficient, and fundamental frequency of the Duffing resonator and the van der Pol oscillator, respectively. 𝑓𝑑 and 𝜔𝑑 are the amplitude and frequency of the external driving force. 𝑘3 is the nonlinearity of the Duffing resonator and 𝑘𝑐 is the coupling stiffness between the two coupled elements. When 𝑘𝑐 diminishes to zero, (1.1) uncouple to yield a driven Duffing resonator and a van der Pol oscillator whose limit cycle is determined by 𝑐2. For the case of 𝑓𝑑=0, (1.1) can be considered as a damped Duffing resonator driven by a van der Pol oscillator. The dynamics of the undriven system has been numerically studied using three control parameters, namely, two damping coefficients (𝑐1 and 𝑐2) and coupling stiffness (𝑘𝑐) [11]. Three different synchronization phenomena were found and a chaotic state was clearly identified in the phase diagram. The dynamics of such an undriven system was also studied in the form of a van der Pol oscillator with a nonlinear restoring force coupled to a Duffing resonator through a mixed velocity and displacement coupling condition [12] and an inertial force condition [13]. Analytic solutions of stable oscillation states were derived in both cases and the chaotic behaviour was also observed. As a result, we can envisage that more complex dynamics could appear in such an externally driven nonlinear system of two elastically coupled oscillators of different types of attractor.

This paper is organized as follows. In Section 2, we give an analytic treatment of (1.1). The method of multiple time scales is used to find approximate solutions of the oscillatory states. We end this section by giving a stability analysis of the fixed points on the stable oscillation curve. Section 3 presents bifurcation analysis of the coupled system in the parameter space. Section 4 is concerned with asymptotic dynamics of the coupled system using a direct integration method. Concluding ideas are given in Section 5.

2. Analytic Treatment

Usually, MEMS resonators and oscillators are working in a high 𝑄 scenario while they are driven into their resonances with a small force. However, if the Duffing resonator and van der Pol oscillator are working at different fundamental frequencies (the nonresonant case 𝜔1𝜔2), it can be proven that the dynamic systems are uncoupled and the time evolution of the amplitudes is the same as the classical driven Duffing resonator and van der Pol oscillator. In the following section, we consider only the case where both the internal and external resonances coincide (𝜔1=𝜔2). As the coupling stiffness varies across a large range, we deal separately with the system of weak coupling (𝑘𝑐1) and strong coupling (𝑘𝑐1).

2.1. The Resonant Case of Weak Coupling

When the van der Pol oscillator is weakly connected to the Duffing resonator, we rewrite (1.1) as follows: ̈𝑥1+𝜖𝜇1̇𝑥1+𝑥1+𝜖𝛼𝑥31𝑥=𝜖𝛽2𝑥1+𝜖𝐹cos(Ω𝜏),̈𝑥2+𝜖𝜇2𝑥221̇𝑥2+𝑥2𝑥=𝜖𝛽1𝑥2,(2.1) where 𝑐1=𝜖𝜇1, 𝑐2=𝜖𝜇2, 𝑘3=𝜖𝛼, 𝑘𝑐=𝜖𝛽, 𝑓𝑑=𝜖𝐹, Ω=𝜔𝑑/𝜔1 and 𝜖 is a small parameter.

The stable response of the system driven near the primary resonance can be calculated from (2.1) using the standard two-time scales method [14]. Briefly, we introduce a fast scale 𝜏0 (𝜏0=𝜏) and a slow scale 𝜏1 (𝜏1=𝜖𝜏0) characterizing the modulation in amplitudes and phases and seek solutions of (2.1) that are expressed in the form of an asymptotic expansion as follows: 𝑥1=𝑥10𝜏0,𝜏1+𝜖𝑥11𝜏0,𝜏1,𝑥2=𝑥20𝜏0,𝜏1+𝜖𝑥21𝜏0,𝜏1.(2.2)

Substituting (2.2) into (2.1) and equating coefficients of like powers of 𝜖, we obtain that𝐷20𝑥10+𝑥10𝐷=0,(2.3a)20𝑥20+𝑥20𝐷=0,(2.3b)20𝑥11+𝑥11=2𝐷0𝐷1𝑥10𝜇1𝐷0𝑥10𝛼𝑥310𝑥+𝛽20𝑥10+𝐹cosΩ𝜏0𝐷,(2.3c)20𝑥21+𝑥21=2𝐷0𝐷1𝑥20𝜇2𝐷0𝑥20𝑥220𝑥1+𝛽10𝑥20,(2.3d)where 𝐷0=𝜕/𝜕𝜏0 and 𝐷1=𝜕/𝜕𝜏1.

The general solution of (2.3a) and (2.3b) can be expressed as 𝑥10=𝐴1𝜏1𝑒𝑖𝜏0𝑥+𝑐.𝑐,20=𝐴2𝜏1𝑒𝑖𝜏0+𝑐.𝑐,(2.4) where 𝐴1 and 𝐴2 are arbitrary constants and 𝑐.𝑐 denotes the complex conjugate.

Substituting 𝑥10 and 𝑥20 into the right-hand side of (2.3c) and (2.3d), we obtain the following set of equations by imposing secular conditions: 2𝐴1+𝜇1𝐴1||𝐴3𝑖𝛼1||2𝐴1𝐴+𝑖𝛽2𝐴1𝑖2𝐹𝑒𝑖𝜎𝜏1=0,2𝐴2+𝜇2||𝐴2||2𝐴12𝐴+𝑖𝛽1𝐴2=0.(2.5) We have introduced the frequency detuning 𝜎 to characterize the closeness of driving frequency to the fundamental frequency, given by Ω=1+𝜖𝜎.

Expressing 𝐴1 and 𝐴2 in polar form: 𝐴1=𝑎1𝑒𝑖𝜃1,𝐴2=𝑎2𝑒𝑖𝜃2,(2.6) and substituting the polar expressions into (2.5) and separating real and imaginary parts, we obtain the following coupled first-order differential equations: 2𝑎1=𝜇1𝑎1+𝛽𝑎2sin𝛾2𝐹2sin𝛾1,2𝑎2=𝜇2𝑎22𝑎12𝛽𝑎1sin𝛾2,2𝑎1𝛾1=(𝛽2𝜎)𝑎1+3𝛼𝑎31𝛽𝑎2cos𝛾2𝐹2cos𝛾1,2𝑎2𝛾2=(𝛽2𝜎)𝑎2𝛽𝑎1cos𝛾22𝑎2𝛾1,(2.7) where 𝛾1=𝜃1𝜎𝜏1 and 𝛾2=𝜃2𝜃1.

To determine the steady-state motion, we use the condition that 𝑎1=𝑎2=0 and 𝛾1=𝛾2=0. From (2.7), we obtain the following set of algebraic equations, 𝑎1(𝛽2𝜎)+3𝑎31𝑎𝛼22(𝛽2𝜎)𝑎12+𝑎1𝜇1𝑎22𝜇2𝑎11𝑎222=14𝐹2,(2.8a)(𝛽2𝜎)2𝑎22+𝜇22𝑎221𝑎222=𝑎21𝛽2.(2.8b)Equation (2.8b) shows Duffing-type behaviour of the van der Pol oscillator, provided that the coupling and damping coefficients are not zero. If the coupling effect diminishes, then (2.8a) will be reduced to the frequency-amplitude modulation of a normal Duffing resonator.

2.2. The Resonant Case of Strong Coupling

If the coupling stiffness 𝑘𝑐 is not the same order of magnitude as the other parameters, the analysis in Section 2.1 does not hold, so a modified analysis considering strong coupling is given as follows. When 𝑘𝑐 is larger than 4.5 for example, 1/(2𝑘𝑐+1) is less than 0.1. Therefore, We choose 𝜖=1/(2𝑘𝑐+1) and assume all other coefficients are small as defined in Section 2.1. Then (1.1) can be rewritten as ̈𝑥1+̈𝑥2+𝜖𝜇1̇𝑥1+𝜖𝜇2𝑥221̇𝑥2+𝑥1+𝑥2+𝜖𝛼𝑥31𝜖=𝜖𝐹cos(Ω𝜏),̈𝑥1̈𝑥2+𝜖2𝜇1̇𝑥1𝜖2𝜇2𝑥221̇𝑥2+𝑥1𝑥2+𝜖2𝛼𝑥31=𝜖2𝐹cos(Ω𝜏).(2.9)

Repeating the same analysis in Section 2.1, we obtain the solutions of (2.9) for the zero-order approximation as 𝐷20𝑥10+𝑥20+𝑥10+𝑥20𝑥=0,(2.10)10𝑥20=0.(2.11) Equation (2.11) means that the Duffing resonator and the van der Pol oscillator are completely synchronised. As shown in Section 2.1, we can easily write down the solutions of (2.10) and (2.11) as follows: 𝑥10=𝑥20=𝐴𝑒𝑖𝜏+𝑐.𝑐.(2.12)

The first-order approximation gives 𝐷20𝑥11+𝑥21+𝑥11+𝑥21=2𝐷1𝐷0𝑥10+𝑥20𝜇1𝐷0𝑥10𝜇2𝑥220𝐷10𝑥20𝛼𝑥310𝑥+𝐹cos(Ω𝜏),11𝑥21=0.(2.13)

Apparently, the Duffing resonator and the van der Pol oscillator are still completely synchronised in the first-order approximation when they are strongly coupled. Substituting the expressions of 𝑥10 and 𝑥20 into the right-hand side of (2.13), we obtain the following equation by imposing secular conditions:4𝐴+𝜇1𝜇2𝐴+𝜇2𝐴||𝐴||2||𝐴||3𝑖𝛼𝐴21+𝑖2𝐹𝑒𝑖𝜎𝜏1=0.(2.14)

Expressing 𝐴 in polar form as defined in Section 2.1 (note that 𝑎1=𝑎2=𝑎, 𝜃1=𝜃2 as (2.12) holds), and substituting into (2.14), we obtain the following set of first-order differential equations for the amplitude and phase: 4𝑎=𝜇2𝜇1𝑎𝜇2𝑎312𝐹sin𝛾1,4𝑎𝛾1=4𝑎𝜎+3𝛼𝑎312𝐹cos𝛾1.(2.15) For the steady state solutions, we obtain the following algebraic equation describing the amplitude response of the system:𝑎24𝜎+3𝛼𝑎22+𝑎2𝜇1𝜇2+𝜇2𝑎22=14𝐹2.(2.16)

After eliminating the secular terms in the two-time scales analysis, the second-order differential equation, (2.13) can be solved. Accordingly, 𝑥11 and 𝑥21 are obtained as follows: 𝑥11=𝑥21=116𝛼𝜇1𝑖𝐴3𝑒3𝑖𝜏+𝑐𝑐.(2.17) Therefore, the general solution of the strongly coupled Duffing resonator and van der Pol oscillator in the first-order approximation is given by𝑥1=𝑥2𝑥=10+𝜖𝑥11=𝑎𝑒𝑖(Ω𝜏+𝛾)+1𝜖16𝛼𝜇2𝑖𝑎3𝑒3𝑖(Ω𝜏+𝛾),+𝑐𝑐(2.18) which shows that in the case of strong coupling, the van der Pol oscillator is completely synchronized to the Duffing resonator and their oscillations are also locked to the external driving force by a phase lag.

2.3. Stability Analysis

The frequency-amplitude relation for weak coupling and strong coupling is, respectively, governed by (2.8a)-(2.8b) and (2.16). By varying the frequency detuning 𝜎, the amplitudes (𝑎1 and 𝑎2) are computed numerically with other parameters provided. Figure 1 shows a frequency response for a weakly coupled system. The amplitude of the Duffing resonator bends towards the higher frequency side like a normal Duffing resonator with a hardening cubic stiffness, but there is another antiresonance peak in addition to the hysteresis domain. Interestingly, even though the van der Pol oscillator is weakly connected to the Duffing resonator, its dynamic behaviour is strongly affected by the Duffing resonator as shown in Figure 1.

The van der Pol response curve reveals a new behaviour. Its resonance is sharpened and a shoulder has been introduced at the location of the antiresonance of the Duffing resonator, but more relevant is the appearance of additional multiple hysteresis branches, denoting new characteristics of the weakly-coupled van der Pol oscillator.

Are these branches stable or not? To answer this, the eigenvalues of the system's Jacobian matrix 𝐽𝑎𝑐 were evaluated at a pair of equilibrium points (𝑎𝑒1,𝑎𝑒2) which correspond to a pair of points in the frequency response curves of Figure 1 for the Duffing resonator and van der Pol oscillator, respectively. The eigenvalues of 𝐽𝑎𝑐 determine linear stability of the equilibrium, which is asymptotically stable if all eigenvalues have negative real parts whereas it is unstable if at least one eigenvalue has positive real part. Figure 2 shows the stable (solid lines) and unstable (dotted lines) branches of the coupled system. As can be expected, the middle branches of the hysteresis for the Duffing resonator and van der Pol oscillator, respectively, are unstable for the system parameters that have been chosen. The system also displays interesting phenomena such as an antiresonance at Ω1 and a shoulder in the curve which is due to a stagnant response of the van der Pol oscillator. The stability of the system is further investigated numerically in the next section where an analysis of the time series is carried out using tools to study chaos in nonlinear dynamical systems. The oscillating states of the system at specific points on the branches are shown in Figure 9.

For the case of strong coupling, the van der Pol oscillator is synchronized to the Duffing resonator and the coupled system oscillates like a single unit, exhibiting Duffing behaviour. This can be illustrated in Figure 3 where the frequency response curve of the system is plotted to first-order approximation.

Similarly, the stability was checked by calculating the eigenvalues of the corresponding Jacobian matrix of (2.15). Choosing several fixed points on one of the branches shown in Figure 3, we found that as in the case of normal Duffing resonators, the top and bottom branches are stable, while the middle branch is unstable. This is further verified by the transient simulation given in Section 4, where transitions between branches are observed.

3. Bifurcation Analysis

In the two-time scales analysis in Section 2, most coefficients of (1.1) are assumed small and of the same order. However, it is very interesting to know the dynamics of the system when the parameters vary in different ranges. In this section, the dependence of the system dynamics on its parameters are studied using a bifurcation analysis. We transform (1.1) into a set of first-order differential equations as follows: ̇𝑥1=𝑦1,̇𝑦1=𝑐1𝑦1𝑥1𝑘3𝑥31+𝑘𝑐𝑥2𝑥1+𝑓𝑑cos(𝑧),̇𝑥2=𝑦2,̇𝑦2=𝑐2𝑥22𝑦12𝑥2+𝑘𝑐𝑥1𝑥2,̇𝑧=Ω,(3.1) where 𝑧=Ω𝜏=(𝜔𝑑/𝜔1)𝜏 and the internal resonance condition 𝜔1=𝜔2 is applied.

Usually, it is easy to adjust the driving force to change the behaviour of the coupled system once it is designed and physically constructed. Therefore, the frequency and the amplitude of the driving force are used as main control parameters in the bifurcation. For better visualisation of the attractors and their bifurcations the dynamics is investigated in the Poincaré section defined by, =𝑥1,𝑦1,𝑥2,𝑦2,𝑧𝑅4×𝑆1𝑧=const..(3.2) In order to investigate the dependence of the system on a single control parameter, several bifurcation diagrams have been computed.

Each bifurcation diagram is calculated by changing the value of the control parameter in small steps. The last computed values (for a particular control parameter value) are used as a new set of initial conditions for the next control parameter value. Equation (3.1) with different driving frequencies is integrated using a fourth-order Runge-Kutta method with the step size of (2𝜋/Ω)/40, where 2𝜋/Ω is the period of the driving force. Numerical solutions corresponding to the first 500 driving cycles are removed as transient. The given bifurcation diagrams show the projection of the attractors in the Poincaré section onto the 𝑥1 and 𝑥2 coordinates versus the control parameter.

Figure 4 shows the symmetry breaking bifurcation of the periodically driven Duffing resonator coupled to a van der Pol oscillator. When the driving frequency Ω varies from 0.6 to 0.85, the system undergoes a period-one to a period-two and then to a period-one behaviour again. If the frequency of the driving force varies from 1.2 to 1.45, the coupled system suddenly enters into a chaotic state after a short range of period-one motion as shown in Figure 5 and then the system transits into period-one motion. There is a narrow window in the chaotic regime showing multiple period motion. If the damping coefficients 𝑐1 and 𝑐2 are decreased, a few more enlarged windows open up as shown in Figure 6, where the system has undergone period-one, chaotic, and multiple-period motion when the frequency varies. The transition between these different dynamics is justified by examining the Lyapunov exponents of the system shown in Figure 7, which shows the weak degree of chaos in specific windows of the driving frequency.

For a single Duffing resonator, there is a critical value of the driving force which can drive it from the simple regime to the chaotic regime. The external force will also influence the van der Pol oscillator through the elastic coupling. As shown in Figure 6, the system experiences a multiple-periods motion at the frequency 1.325. Using the same parameters, we investigated the effect of the driving force by varying its amplitude from 0.1 to 10. As shown in Figure 8, the system exhibits complex behaviour when the force changes.

4. Dynamics of the Coupled System

The dynamics of the coupled system was studied by integrating (1.1) directly using the fourth-order Runge-Kutta method. As indicated in Section 2.1, the weakly coupled Duffing-van der Pol system has complicated dynamics. Using the parameters in Figure 1 and choosing the frequency Ω at 1.01, where stable and unstable branches of the frequency response curves coexist, we computed the transient motion of the Duffing resonator and van der Pol oscillator. The Poincaré section in Figure 9 shows the system converges to an attractor after a long transient motion and the system shows weak chaos with the largest Lyapunov exponent being about 0.007.

When the coupling stiffness becomes stronger, the van der Pol oscillator quickly synchronizes to, and behaves like the Duffing resonator. The numerical experiment of sweeping the frequency of the driving force clearly illustrates a jumping/hysteresis phenomenon as shown in Figure 10, which is further confirmation of the two-times-scale analysis in Section 2.

5. Conclusion

In this paper, we have studied the dynamics of a coupled nonlinear system consisting of a van der Pol oscillator coupled to a periodically driven Duffing resonator, both of which show different attractors individually. The competition between the two different attractors depends upon the coupling constant and the control parameters.

In the strong coupling regime, dominating synchronization is observed, where the van der Pol oscillator is locked to the behaviour of the Duffing resonator.

The weak coupling regime reveals a multitude of complex behaviours, including multiperiod transitions and weakly chaotic motion which was confirmed by an analysis of Lyapunov exponents of the system Jacobian. We have also observed, multistable branches in the asymptotic amplitude-frequency response curves, and also an antiresonance in the Duffing response leading to a sharpening and a shoulder in the van der Pol response curve.

These unexpected characteristics will have profound consequences when we construct chains and arrays composed of repeated units of coupled Duffing and van der Pol systems. Arrays of such components could have practical use as compact and resilient sensor arrays, for example. However the response of a sensor array of such coupled subsystems is likely to have a complexity of response characteristics beyond normal expectations of current sensor arrays. The next challenge will be in designing sensor arrays based on these subsystems which utilise this vast spectrum of computational behaviours.


The authors acknowledge the financial support of this work from the Engineering and Physical Sciences Research Council (EPSRC) under the Project nos. EP/D00036X/1 and EP/D501032/1.