Stabilization of dynamics of oscillatory systems by nonautonomous perturbation

Synchronization and stability under periodic oscillatory driving are well understood, but little is known about the effects of aperiodic driving, despite its abundance in nature. Here, we consider oscillators subject to driving with slowly varying frequency, and investigate both short-term and long-term stability properties. For a phase oscillator, we find that, counterintuitively, such variation is guaranteed to enlarge the Arnold tongue in parameter space. Using analytical and numerical methods that provide information on time-variable dynamical properties, we find that the growth of the Arnold tongue is specifically due to the growth of a region of intermittent synchronization where trajectories alternate between short-term stability and short-term neutral stability, giving rise to stability on average. We also present examples of higher-dimensional nonlinear oscillators where a similar stabilization phenomenon is numerically observed. Our findings help support the case that in general, deterministic nonautonomous perturbation is a very good candidate for stabilizing complex dynamics.

Figure 1

(a) Time-varying existence of the attracting point of Eq. (5). The $\left[\mathrm{\Delta }\omega \right(t\left)\right]$-dependent curve of $\dot{\psi }$ against $\psi$ moves up and down over time, as indicated by the two solid lines representing where the curve could be at two different instants in time. When this curve lies between the dashed lines, the system has an attracting point, and otherwise, not. (b) Phase diagram showing three regimes. Light (region I in the text), medium (region III), and dark gray (region II) show where the system is never synchronizing, intermittently synchronizing, and always synchronizing, respectively. Solid white curves show the border between synchronization and nonsynchronization for if $f(·)$ is set to 0; and white dashed curves show the border between synchronization and nonsynchronization for if $f(·)$ is set to $\pm 1$. When $k$ is increased, regions I and II decrease while region III increases; hence, in particular, the Arnold tongue consisting of the union of regions II and III increases.

Figure 2

Numerically obtained long-term maximum Lyapunov exponent ${\lambda }_1$ over parameter space for (3), with different $k$. The LE are computed over five cycles of the frequency modulation (about 630 s). In each case, 20 random initial conditions were taken from the square $[-1,1]\times [-1,1]$, and the average maximum LE over these trajectories is plotted. (a) $k=0$, (b) $k=0.1$, (c) $k=0.4$, (d) $k=0.8$. The Arnold tongue (shades of blue) is enlarged as $k$ increases. Gray represents zero values.

Figure 3

Trajectories of five solutions ${\theta }^1\left(t\right),...,{\theta }^5\left(t\right)$ of Eq. (1), with initial conditions ${\theta }^i\left(0\right)=\frac{i-1}{5}2\pi$, subject to the same driving ${\theta }_1\left(t\right)$. Here, $\gamma =2.5$ and $k=0.4$. In (a), ${\omega }_0=4$, and so the system is in region II according to Eq. (9); in (b), ${\omega }_0=6$, and so the system is in region III according to Eq. (10); in (c), ${\omega }_0=9$, and so the system is in region I according to Eq. (8). In each case, the upper plot shows the first 8 s of the sine of the five trajectories, while the lower plot shows the distance between ${\theta }^1\left(t\right)$ and ${\theta }^2\left(t\right)$ over about the first 630 s (more precisely, five cycles of the frequency modulation); in (a) and (b), the inner graph shows the same information on a logarithmic scale. In (a) and (b), the system lies within the Arnold tongue as described in Sec. 2c, and the five trajectories are observed to synchronize and to remain in synchrony; by contrast, in (c), the system does not lie within the Arnold tongue, and no synchronization is observed.

Figure 4

Analysis of the time-variable dynamical properties of a trajectory of (3), in the three regions (see Fig. 1); in (a)–(f), a random initial condition is taken from the square $[-1,1]\times [-1,1]$. (a)–(c) Time-frequency representation (showing amplitude) of $\theta \left(t\right)$ obtained as the (unwrapped) polar angle of the trajectory of (3), extracted using a continuous Morlet wavelet transform ($p=1$) with central frequency 3; and (d)–(f) maximum FTLE for the trajectory of (3) over a time window of width $\tau =0.1$ s for regions II, III, and I, from left to right. Parameters are set to $k=0.4$ and $\gamma =2.5$; in (a) and (d) ${\omega }_0=4$, in (b), (e), and (g) ${\omega }_0=6$, and in (c), (f), and (h) ${\omega }_0=9$. (a), (d) Region II exhibits frequency entrainment and a stable phase at all times. (b), (e) Region III shows intermittent, but regular, epochs of frequency entrainment; the phase is stable on average over a long time. (c), (f) Region I, no frequency entrainment, and a FTLE rapidly oscillating around zero. (g), (h) Prediction (red) of the main observed frequency of the system over time, based on Eq. (12) with values of ${\mathrm{\Omega }}_{\psi }\left(t\right)$ taken from the ${\mathrm{\Omega }}_{\psi }$ curve with $k=0$ in Fig. 6, for (g) region III and (h) region I. Interestingly, in region I, the main observed frequency oscillates in antiphase with those of the driving frequency (dashed). Note that we here only predict the main frequency, and not the higher harmonics observed in (b) and (c).

Figure 5

Analytical and numerical characterization of the three regions based on the proportion of time spent in the stable regime, over $(\gamma ,{\omega }_0)$-parameter space with $k=0.4$. In (a), and (c), $P_t$ is calculated numerically based on time-localized maximum LE (with window $\tau =0.1$ s) for 20 trajectories of (3) with random initial conditions in the square $[-1,1]\times [-1,1]$, over four cycles of the frequency modulation (about 500 s), and the result is averaged over the 20 trajectories. In (b) and (d), the analytical result according to Eq. (11) is shown. Plots (c) and (d) show $P_t$ from (a) and (b) for selected $\gamma$ values. Three distinct regions appear clearly, respectively with $P_t$ values of 0 (region I), 1 (region II), and in between 0 and 1 (region III). Analytical and numerical characterizations show good agreement.

Figure 6

Numerically obtained time-averaged stability properties for a trajectory of Eq. (1) starting at $\theta \left(0\right)=0$, computed over ten cycles of the frequency modulation (about 1260 s). (a) Lyapunov exponent ${\lambda }_1$, and (b) average frequency difference ${\mathrm{\Omega }}_{\psi }$, for $\gamma =2.5$ and different values of the frequency modulation amplitude $k$. For $k=0$, the region of phase stability (as given by ${\lambda }_1<0$) and the region of permanent frequency entrainment (as given by ${\mathrm{\Omega }}_{\psi }=0$) coincide; but for $k>0$, the regions do not coincide: as $k$ is increased, the region with ${\lambda }_1<0$ is increased while the region with ${\mathrm{\Omega }}_{\psi }=0$ is decreased. The graphs look very similar to those in the case of harmonic driving with bounded noise [5]; therefore, investigation of time-variable dynamics for nonautonomously driven oscillators is necessary for an accurate understanding of the dynamical nature of the system.

Figure 7

Numerically obtained time-averaged stability properties for autonomous (plain black), nonautonomous (dashed red), and bounded noise (dot-dashed black) driving, computed over about 1260 s (ten cycles of the periodic frequency modulation used for the nonautonomous case). For the autonomous and nonautonomous cases, Eqs. (4) and (5), respectively, were integrated, with $\psi \left(0\right)=0$; for the noisy case, one sample path of $\xi \left(t\right)$ was generated, and Eq. (13) was integrated with $\psi \left(0\right)=0$, using the same sample path $\xi \left(t\right)$ for all ${\omega }_0$ values. Parameters are set to ${\omega }_1=4$ and $\gamma =1$; for the nonautonomous case, $f=sin(·), k=0.1$, and ${\omega }_m=0.05$, and for the noisy case $D=1.6$ and $\mu =10$. (a) Lyapunov exponent ${\lambda }_1$ and (b) average frequency difference ${\mathrm{\Omega }}_{\psi }$. The nonautonomous and noisy cases, observed on average, present the same enlarging of the negative Lyapunov exponent region, and their ${\mathrm{\Omega }}_{\psi }$ is almost exactly identical, including the plateau.

Figure 8

Dynamics of driven oscillator $\theta$ with (a)–(d) nonautonomous and (b), (e), (f) bounded noise driving; here $\theta$ evolves according to Eq. (1) with ${\theta }_1\left(t\right)={\omega }_1[t-\frac{k}{{\omega }_m}cos\left({\omega }_{m}t\right)]$ in the nonautonomous case, ${\theta }_1\left(t\right)={\omega }_{1}t-{\int }_0^t\xi \left(s\right)ds$ in the noisy case, and $\theta \left(0\right)={\theta }_1\left(0\right)$ in both cases. Parameters are set to ${\omega }_1=4, \gamma =1, {\omega }_0=3$; for the nonautonomous case $k=0.1$ and ${\omega }_m=0.05$, and for the noisy case $D=1, \mu =4$. First $\psi \left(t\right)$ is numerically obtained by integrating Eq. (5) or (13) as appropriate, with $\psi \left(0\right)=0$, and then $\theta \left(t\right)$ is obtained by $\theta \left(t\right)=\psi \left(t\right)+{\theta }_1\left(t\right)$. (a) Sine of the driven phase $\theta$, in the nonautonomous setting. (b) Phase difference $\psi$ between the driving and driven oscillators, over time. The nonautonomous case presents epochs of phase locking as seen by the regular plateaus, whereas the noisy phases' difference drifts without ever phase locking to the driving, with an average velocity that is close to that of the nonautonomous curve. (c), (e) Time-frequency representation (showing power, i.e., square of the amplitude) for $\theta$ extracted using continuous Morlet wavelet transform (with $p=1$) with central frequency 3, and (d), (e) the associated time-averaged power. The main difference is the presence of intermittent frequency entrainment in the nonautonomous case. The average-power spectra are very similar, and do not clearly distinguish the two cases. Long-term LEs were also found to be negative in both cases: for the nonautonomous Eq. (5) with $\psi \left(0\right)=0$, the LE over 500 s was about $-0.32$, and for the noisy Eq. (13) with $\psi \left(0\right)=0$, the LE over 500 s was about $-0.24$.

Figure 9

Numerically obtained long-term maximum Lyapunov exponent ${\lambda }_1$ over parameter space for (3) with ${\theta }_1\left(t\right)={\omega }_1\{t+\frac{0.5k}{{\omega }_m}[sin\left({\omega }_{m}t\right)+\frac{4}{\pi }sin({\omega }_m\pi t/4)]\}$. Here, as in Sec. 2d, ${\omega }_1=4, {\omega }_m=0.05, r_p=1$, and $\epsilon =5$. The LE are computed over five cycles of the frequency modulation (about 630 s). In each case, 20 random initial conditions were taken from the square $[-1,1]\times [-1,1]$, and the average maximum LE over these trajectories is plotted. (a) $k=0$, (b) $k=0.4$, (c) $k=0.8$. The Arnold tongue (shades of blue) is enlarged as $k$ increases. Gray represents zero values.

Figure 10

Numerically obtained long-term maximum Lyapunov exponent ${\lambda }_1$ over parameter space for the forced weakly nonlinear vdP oscillator [see Eq. (15)], with $\epsilon =0.1, {\omega }_1=1$, and ${\omega }_m=0.02$, for different amplitude of frequency modulation $k$. In each case, 20 random initial conditions $[x\left(0\right),\dot{x}\left(0\right)]$ were taken from the square $[-1,1]\times [-1,1]$, and the average maximum LE over these trajectories is plotted. The negative LE region (blue shades) increases as $k$ is increased. Gray represents zero values.

Figure 11

Intermittency in the typically forced vdP system (15). Parameters are set to $k=0.5, \gamma =1, {\omega }_0=1, \epsilon =0.1, {\omega }_m=0.02, {\omega }_1=1$; results are for a random initial condition $[x\left(0\right),\dot{x}\left(0\right)]$ taken from the square $[-1,1]\times [-1,1]$. (a) Time-frequency representation (showing amplitude) of $x\left(t\right)$, extracted using a continuous Morlet wavelet transform ($p=1$) with central frequency 2. (b) Shorter-time-window FTLE max, with window length 0.1 s, is shown in black, and longer-time-window FTLE max, with window length 25 s, is shown in red. The longer-time-window FTLE alternates between epochs where it is negative, and epochs where it is zero. These epochs of negative values coincide with those epochs where, in (a), there appears to be a single main peak in the instantaneous power-frequency spectrum.

Figure 12

Numerically obtained long-term maximum Lyapunov exponent ${\lambda }_1$ over parameter space for the diffusively phase forced van der Pol oscillator [see Eq. (16)], with $\epsilon =0.1, {\omega }_1=1$, and ${\omega }_m=0.02$, for different amplitude of frequency modulation $k$. In each case, 20 random initial conditions $[x\left(0\right),\dot{x}\left(0\right)]$ were taken from the square $[-1,1]\times [-1,1]$, and the average maximum LE over these trajectories is plotted. The stability region (shades of blue) is enlarged as $k$ increases, and chaotic (red) points in (a) are turned stable (shades of blue) in (c). Gray represents zero values.

Figure 13

Numerically obtained long-term maximum Lyapunov exponent ${\lambda }_1$ over parameter space for the coupled and forced Duffing oscillator (17), for different amplitude of frequency modulation $k$. In each case, 20 random initial conditions $[x\left(0\right),\dot{x}\left(0\right),x_1\left(0\right),{\dot{x}}_1\left(0\right)]$ were taken from ${[-1,1]}^4$, and the average maximum LE over these trajectories is plotted. The region of chaotic behavior (red) is reduced as many points are turned stable (shades of blue) as $k$ is increased. Gray represents zero values.