Non-feedback technique to directly control multistability in nonlinear oscillators by dual-frequency driving

Non-feedback technique to directly control multistability in nonlinear oscillators by... Nonlinear Dyn https://doi.org/10.1007/s11071-018-4358-z ORIGINAL PAPER Non-feedback technique to directly control multistability in nonlinear oscillators by dual-frequency driving GPU accelerated topological analysis of a bubble in water Ferenc Hegedus ˝ · Werner Lauterborn · Ulrich Parlitz · Robert Mettin Received: 28 December 2017 / Accepted: 12 May 2018 © The Author(s) 2018 Abstract A novel method to control multistability ing the simulations, the control parameters are the two of nonlinear oscillators by applying dual-frequency amplitudes of the acoustic driving at fixed, commensu- driving is presented. The test model is the Keller– rate frequency pairs. The high-resolution bi-parametric Miksis equation describing the oscillation of a bub- scans in the control parameter plane show that a period- ble in a liquid. It is solved by an in-house initial-value 2 attractor can be continuously transformed into a problem solver capable to exploit the high computa- period-3 one (and vice versa) by proper selection of tional resources of professional graphics cards. Dur- the frequency combination and by proper tuning of the driving amplitudes. This phenomenon has opened Electronic supplementary material The online version of a new way to drive the system to a desired, pre-selected this article (https://doi.org/10.1007/s11071-018-4358-z) attractor directly via a non-feedback control technique contains supplementary material, which is available to without the need of the annihilation of other attractors. authorized users. Moreover, the residence in transient chaotic regimes can also be avoided. The results are supplemented with F. Hegedus ˝ Department of Hydrodynamic Systems, Faculty of simulations obtained by the boundary-value problem Mechanical Engineering, Budapest University of solver AUTO, which is capable to compute periodic Technology and Economics, Budapest, Hungary orbits directly regardless of their stability, and trace e-mail: fhegedus@hds.bme.hu them as a function of a control parameter with the W. Lauterborn · R. Mettin pseudo-arclength continuation technique. Cavitation Bubble Dynamics Group, Drittes Physikalisches Institut, Georg-August-Universität Göttingen, Göttingen, Keywords Control of multistability · Dual-frequency Germany e-mail: werner.lauterborn@phys.uni-goettingen.de driving · Bifurcation structure · GPU programming · Keller–Miksis equation · Nonlinear dynamics R. Mettin e-mail: robert.mettin@phys.uni-goettingen.de U. Parlitz ( ) 1 Introduction Research Group Biomedical Physics, Max Planck Institute for Dynamics and Self-Organization, Göttingen, Germany e-mail: ulrich.parlitz@ds.mpg.de One of the common features of nonlinear systems is multistability [1]; that is, two or more attractors can U. Parlitz coexist at a fixed parameter combination. Multistabil- Institute for Nonlinear Dynamics, ity appears in almost any field of science; for instance in Georg-August-Universität Göttingen, Göttingen, Germany biological systems [2,3], hydrodynamics [4], mechan- 123 F. Hegedus ˝ et al. ical engineering [5], chemical reactions [6,7], neuron There are three main control strategies known in the dynamics [8], climate dynamics [9] or social systems literature, namely non-feedback, feedback and stochas- [10] to name a few. Multistability is also found in bub- tic control [1]. Each has its own advantages and disad- ble dynamics and acoustic cavitation [11], thetopicof vantages. The non-feedback techniques are simple and the present study. In applications, a major challenge in easy to use, where the main aim is to kick the system nonlinear dynamics is how to control multistability and to another stable state [31] or annihilate attractors by switch between coexisting states representing different periodic perturbation (modulation) of a parameter or system behaviour. a state variable [32]. Unfortunately, it cannot be guar- The mechanism of the emergence of multiple attrac- anteed that the system settles down onto the desired tors can be rather different. In conservative systems, for attractor. The kick introduces randomness to the con- example, one can obtain an arbitrarily large number of trol, and extremely long transients are also possible coexisting attractors by introducing small dissipation in the presence of transient chaos. In case of periodic [12] (see also the Kolmogorov–Arnold–Moser (KAM) modulation, there is no direct control over which attrac- theory [13]). The number of attractors is roughly tors are annihilated, since usually the attractor with the inversely proportional to the damping factor [14]. The smallest basin is destructed first. With increasing per- appearance of homoclinic tangencies or the coupling turbation magnitude, the annihilation process contin- of systems with a large number of unstable invariant ues with attractors having successively the next smaller manifolds can also lead to a high level of multistabil- basin. Applying a feedback control strategy, the desired ity. In these cases, even an infinite number of attracting attractor can be targeted directly; moreover, it is a reli- periodic orbits called Gavrilov–Shilnikov–Newhouse able technique to control attractors against external sinks can coexist [15,16]. Another important mecha- noise. However, details are necessary about the state nism leading to multistability is the coupling of a large space [33] or even the solution of the attractor itself number of identical systems, where the number of the [34] or its Jacobian. This makes the feedback control attractors scales with the number of oscillators [17]. technique unusable in certain problems. For instance, in This results in a so-called attractor crowding that makes an acoustically trapped bubble, it is hard to obtain state the system extremely sensitive to external noise lead- space information about the oscillating bubble and feed ing to random hopping between many coexisting states. this information back to the resonator. The third kind Several other mechanisms such as delayed feedback of control (stochastic) also performs attractor annihila- [18], parametric forcing [19] or external noise [20] can tion, this time with the addition of external noise [35]. also induce multistability. The drawback is similar as in case of the non-feedback The importance of the control of multistability can technique: the attractors are destructed in the order of be explained with the fact that different stable states the size of their basin of attraction [1]. That is, the direct represent different system performances. A good exam- control is again lost. Attractors with very small basins ple is the large subject of chaos control [21–25], where are usually buried by the noise immediately, without the main aim is to eliminate unpredictable behaviour. the need of increasing the magnitude of the noise to Furthermore, there are various fields where the con- some higher level. trol of multistability is of paramount importance: laser The main aim of the present study is to propose a pos- technology [26], optical communication [27], cardiol- sible non-feedback control technique capable of driv- ogy [28], genetics [29] or ecology [30]. If multistability ing a sinusoidally excited system to a pre-selected peri- is undesired, then basically two approaches are possi- odic attractor directly. It is based on the observation that ble: make the system monostable or preselect and sta- a period-2 and a period-3 attractor can be continuously bilize an attractor against external noise. In case of a transformed into each other by adding a second, sinu- desirable multistability, where the control between dif- soidal component to the excitation. The proper tuning ferent states (operations) on demand is important, the of the excitation amplitudes is necessary; moreover, task is to properly select an attractor the system should the fixed commensurate frequency combination must approach and/or exclude certain undesired attractors be adjusted according to the periods of the attractors from the dynamics. In these cases, the stabilization being transformed. The advantages and disadvantages against external noise can also be crucial. over other, already existing control techniques, the lim- itations and possible extensions to transformation of 123 Nonlinear oscillators by dual-frequency driving a period-n to a period-m attractor in general, and the and the vapour, p = 3166.8Pa. p denotes the pres- V L application to other nonlinear oscillators are also dis- sure in the liquid at the bubble-liquid interface. The cussed. surface tension is σ = 0.072 N/m, and the liquid kine- −4 The test model is the Keller–Miksis equation well matic viscosity is μ = 8.902 Pa s. The gas inside known in the bubble dynamics community. It is a the bubble obeys a simple polytropic relationship second-order nonlinear ordinary differential equation 3γ describing the oscillation of a bubble in a liquid [11]. 2σ R p = P − p + , (4) G ∞ V The parameter space of the present study is relatively R R large due to the dual-frequency driving (amplitudes and frequencies). The approach to obtain a global where the polytropic exponent is γ = 1.4 (adiabatic behaviour), the equilibrium bubble radius is R = picture about the bifurcation structure in the four- dimensional driving parameter space within reason- 10 μm and the static pressure is P = 1 bar. able time is to exploit the high arithmetic processing For numerical purposes, system (1)–(4) is rewritten power of professional graphics cards (GPUs) using an into a dimensionless form by introducing the dimen- in-house initial-value problem (IVP) solver written in sionless variables C++/CUDA C. The parameter scans are supplemented with results obtained by the boundary-value problem τ = t, (5) 2π solver AUTO [36], which is suitable to compute peri- odic orbits directly regardless of their stability, and fol- y = , (6) low their paths with the pseudo-arclength continuation 2π technique. y = R , (7) R ω E 1 2 Mathematical model and the equations are rearranged to minimize the num- ber of coefficients. The final form of the numerical The test model studied (Keller–Miksis equation [11]) model is describing the evolution of the bubble radius reads ˙ ˙ R R 3 ¨ ˙ 1 − R R + 1 − R y ˙ = y , (8) 1 2 c 3c 2 L L KM R R d ( p − p (t )) y ˙ = , (9) L ∞ 2 = 1 + + , (1) D KM c c dt ρ L L L where R(t ) is the time dependent bubble radius; c = where the numerator, N , and the denominator, D , L KM KM 1497.3m/s and ρ = 997.1kg/m are the sound speed are and density of the liquid domain, respectively. Accord- N = (C + C y ) − C (1 + C y ) ing to the general, dual-frequency treatment, the pres- KM 0 1 2 2 9 2 sure far away from the bubble 1 y y 3 2 2 −C − C − 1 − C y 3 4 9 y y 3 2 1 1 p (t ) = P + P sin(ω t ) + P sin(ω t + θ) (2) ∞ ∞ A1 1 A2 2 − (C sin(2πτ ) + C sin(2π C τ + C )) 5 6 11 12 1 + C y − y C cos(2πτ ) consists of a static component, P , and periodic com- ( ) ( 9 2 1 7 ponents with pressure amplitudes P and P , angular A1 A2 + C cos(2π C τ + C )) , (10) 8 11 12 frequencies ω and ω , and with a phase shift θ.The 1 2 and connection between the pressures inside and outside the bubble at its interface can be written as D = y − C y y + C C . (11) KM 1 9 1 2 4 9 2σ R p + p = p + + 4μ , (3) G V L L The coefficients are summarized as follows: R R where the total pressure inside the bubble is the sum of 1 2σ 2π C = P − p + , (12) 0 ∞ V the partial pressures of the non-condensable gas, p , ρ R R ω L E E 1 123 F. Hegedus ˝ et al. 1 − 3γ 2σ 2π tions forward in time and, after the convergence of the C = P − p + , (13) 1 ∞ V ρ c R R ω transient trajectory, record and save the important prop- L L E E 1 erties of the solutions found. Thereby, our main inter- P − p 2π ∞ V C = , (14) est is in periodic solutions and their organization in the ρ R ω L E 1 reduced driving parameter space discussed below. The 2σ 2π first part (2048 T ) of the integration is regarded as ini- C = , (15) ρ R R ω L E E 1 tial transient and dropped followed by integration of 4μ 2π the system further by 8192 T for the proper conver- C = , (16) ρ R ω 1 gence of average quantities like Lyapunov exponent or winding number, where T is the period of the dual- P 2π A1 C = , (17) frequency excitation function. The scheme employed ρ R ω L E 1 is the explicit and adaptive Runge–Kutta–Cash–Karp P 2π A2 method with embedded error estimation of orders 4 C = , (18) ρ R ω L E 1 and 5 (the algorithm is adapted from [38]). The quan- tities stored are the points of the (global) Poincaré sec- ω P 2π 1 A1 C = R , (19) 7 E tion (see the choice in Sect. 4 below), which are stan- ρ c R ω L L E 1 dard ingredients for bifurcation diagrams; furthermore, ω P 2π 1 A2 the maximum bubble radius and the maximum absolute C = R , (20) 8 E ρ c R ω L L E 1 value of the bubble wall velocity, which are important R ω E 1 for applications; finally, the period,the Lyapunov expo- C = , (21) 2π c L nent and the winding number of the attractors found, C = 3γ, (22) quantities that are essential for a detailed analysis of bifurcation structures. C = , (23) A strategy to represent the results of parametric stud- ies involving high-dimensional parameter spaces con- C = θ. (24) sists in creating high-resolution bi-parametric plots, a rapidly spreading technique in the investigation of The angular frequencies ω and ω are normalized 1 2 nonlinear systems with a high-dimensional parame- by the linear, undamped angular eigenfrequency [37] ter space [39–49]. The system studied here, a bub- ble in water with dual-frequency acoustic excita- 3γ( P − p ) 2(3γ − 1)σ ∞ V ω = − (25) tion, has a four-dimensional driving parameter space 2 3 ρ R ρ R L L E E ( P , P ,ω ,ω ). For simplicity, the phase shift is A1 A2 R1 R2 set to θ = 0. It is restricted in the present study by the of the unexcited system that defines the relative fre- requirement that the relative-frequency pairs obtained quencies as from ω and ω are rational. The pressure amplitudes R1 R2 1 P and P are taken as the main control parameters at A1 A2 ω = , (26) R1 ω fixed relative-frequency pairs. The range of each pres- sure amplitude is 0–5 bar, investigated at first with a ω = . (27) R2 resolution of 501 steps for an overview. The relative- frequency values ω and ω are chosen from the R1 R2 With R = 10 µm and with the other constants in (25) following set: given before in this section, the eigenfrequency is f = ω /2π = 340 kHz. 1 1 1 1 1 2 3 5 10 , , , , , , , , . (28) 10 5 3 2 1 1 1 1 1 3 Numerical implementation and parameter choice Exploiting symmetries in the driving parameter space, this gives i = 36 combinations of two fre- i =1 The simplest technique to solve the system (8)–(9)isto quencies. Observe that two orders of magnitude dif- take an initial-value problem solver, integrate the equa- ference in the frequency range are covered and that, 123 Nonlinear oscillators by dual-frequency driving for instance, ω = 1/5 and ω = 1/1 means that R1 R2 T T 1 2 =2 R2 in an experiment, the bubble is driven by 0.2 ω and by the resonance frequency ω . Computations are per- formed at every possible frequency combination of the above set, which means altogether 36 two-dimensional =3 -1 plots for every quantity to be studied (period, Lya- R1 dual-frequency excitation punov exponent, etc). Moreover, owing to the fine -2 0 0.5 1 1.5 2 2.5 3 details in the diagrams, high-resolution surveys are [-] necessary. In order to reveal possible coexisting attrac- tors, 10 randomized initial conditions are used at each Fig. 1 Example for graphs of a dual-frequency excitation (red) and its two harmonic components as mainly used in the study parameter set. Indeed, the 10 different trajectories can (ω = 3, black; ω = 2, blue) with pressure amplitudes R1 R2 converge to different attractors with distinct proper- P = P = 1 bar. The vertical lines marked by T (black A1 A2 1 ties (e.g. period). Altogether, a single plot contains solid), T (blue solid) and T (red dashed) designate the periods 501 × 501 × 10 ≈ 2.5 million initial-value prob- of the high-frequency component, the low-frequency component lems. and the dual-frequency driving, respectively. T = 3 × T = 2 × T = 3. (Color figure online) Even a single bi-parametric plot requires large com- putational capacities; thus, the high processing power of professional videocards (GPGPUs) is exploited. During the parameter studies, millions of simple 4 Global Poincaré section and the period-reducing (only second order) and independent ordinary dif- phenomenon ferential equations are to be solved. Moreover, the algorithm employed needs only function evaluations. The parameters investigated in the present study are Therefore, our problem is well suited for paralleliza- the pressure amplitudes ( P , P ) and the relative fre- A1 A2 tion on GPUs. The program code is written in a quencies (ω ,ω ) of the dual-frequency excitation R1 R2 C++ and CUDA C software environment. The avail- with phase shift θ = 0, i.e. C = 0. The main con- able GPUs are a Titan Black card (Kepler archi- trol parameters are the pressure amplitudes at several tecture, 1707 GFLOPS double-precision processing fixed frequency combinations. The ratio of each fre- power), two Tesla K20 cards (Kepler architecture, 1175 quency pair is rational; therefore, periodicity of the GFLOPS) and eight Tesla M2050 cards (Fermi archi- dual-frequency driving is guaranteed (quasiperiodic tecture, 515 GFLOPS). The application of double- forcing is excluded). precision floating-point arithmetic is mandatory in bub- As an example, a dual-frequency excitation is pre- ble dynamics due to the possibly large-amplitude, sented in Fig. 1 as a function of time together with its collapse-like response of the system [11]. The final, two periodic components at ω = 3 and ω = 2. For R1 R2 optimized code is approximately about 50 times faster definiteness, the pressure amplitudes are set to P = A1 on the Titan Black, about 30 times faster on the Tesla P = 1 bar. The black and blue curves are the high- A2 M2050 and about 120 times faster on the Tesla K20 and low-frequency components, respectively. Since the than on a four-core Intel Core i7-4790 CPU. Sur- dimensionless time τ is defined by means of the first fre- prisingly, the professional Tesla K20 card was more quency component, see Eq. (5), the normalized period than two times faster than the gamer Titan Black card corresponding to the relative frequency ω is T = 1, R1 1 even though the theoretical peak processing power is indicated by the vertical black line in Fig. 1. The nor- higher for the Titan Black GPU. A thorough perfor- malized period of the second component (relative fre- mance analysis between GPUs and CPUs solving large quency ω ) can be calculated with the frequency ratio; R2 numbers of ordinary differential equations is beyond that is, T = 1/C = ω /ω = 3/2 = 1.5 (blue 2 11 R1 R2 the scope of the present study. The interested reader vertical line). Compare these periods with the argu- is referred to the publications [38,50–52]for more ments of the harmonic functions in equation (10) and details. keep in mind that the phase shift is zero (θ = C = 0). As only normalized periods will be used throughout the article, this very often occurring quantity will be simply called “period”. p , p [bar] A1 A2 F. Hegedus ˝ et al. The sum of the two harmonic components results is no restriction to small-amplitude perturbations cor- in a still periodic, but non-harmonic forcing function responding to the second sinusoidal component. presented by the red curve in Fig. 1. The period of this function is T = 3T = 2T = 3, whichisusedasa 5.1 The branching phenomenon 1 2 global Poincaré section for the system (8)–(9). Con- sequently, during the numerical integration of the sys- Out of the data sets for the 36 relative-frequency com- tem with an initial-value problem solver, the continu- binations, four diagram pairs are shown in Fig. 2 at the ous trajectories are sampled at time instants τ = n · T four different relative-frequency pairs (ω ,ω ) = n R1 R2 (n = 0, 1, 2,...). Accordingly, the period, T ,ofa (0.2, 0.1), (1, 0.1), (3, 1) and (3, 2) organized in rows periodic orbit is defined as T = N · T , N ∈ IN, and and for two different quantities (Lyapunov exponent the solution is called period- N orbit. Observe, however, and period) organized in columns. The left column that in the special cases of P = 0 bar or P = 0 bar, contains the Lyapunov-exponent diagrams, where the A1 A2 the period T , which defines the global Poincaré section, colour-coded area means chaotic solutions (positive does not coincide with the period of the actual driving exponent) and where in the greyscale domains there T or T .For P = 0 bar, the dual-frequency forc- are periodic attractors (negative exponent). In the right 2 1 A1 ing becomes a purely harmonic function shown by the column, the periods of the (converged) solutions are blue curve in Fig. 1, and the continuous trajectories presented up to period-6. Periodic solutions higher than are sampled only after every second period of excita- six including chaos can be found in the black domains. tion, 2T . Therefore, an originally period-2k solution In the upper two rows, the plots correspond to relative is observed only as a period-k orbit. Similarly, when frequencies lower than or equal to the linear resonance P = 0 bar, the trajectories are sampled only after frequency (ω ≤ 1), while the lower two rows have A2 R1,2 every third period of the excitation, 3T (compare the relative frequencies higher than or equal to the linear black and red vertical lines in Fig. 1), and the period resonance frequency (ω ≥ 1). R1,2 of an originally period-3k orbit is reduced from 3k to Since detailed investigations of dual-frequency driven k.This period-reducing phenomenon, occurring when systems for arbitrary driving amplitudes with ratio- one of the two components of the forcing function is nal frequency ratio are absent in the literature, only vanishing and only a purely harmonic excitation is left, very few information exists about the bifurcation struc- plays an important role in the bifurcation structure of ture shown in Fig. 2. It can be clearly seen that the the system, discussed in more detail in Sect. 5; further- dual-frequency driving causes a very complex interplay more, period reduction necessarily takes place for other between the resonances originating from the single- frequency pairs as long as the ratio of a pair is rational. frequency driven system (horizontal and vertical axes in Fig. 2), see especially the first two rows of the figure. 5 Results from a global scan Their comprehensive topological analysis is beyond the scope of the present study; instead, it intends to add Throughout the following subsections, the fundamental some pieces to our knowledge by looking for special mechanism capable of controlling multistability in har- features in the diagrams. In particular, we focus on a monically driven nonlinear oscillators is presented. The very specific phenomenon, the branching mechanism, basis is the period-reducing phenomenon discussed turning out to be a special feature of dual-frequency above, by which the dual-frequency treatment can gen- driven nonlinear oscillators with rational frequency erate “bridges” between different periodic attractors. combinations, which can be used for the control of mul- First, a global picture is shown via a series of high- tistability. resolution bi-parametric plots of Lyapunov-exponent Observe that on the horizontal axis in the third row and periodicity diagrams at different frequency com- and second column of Fig. 2, three period-3 branches binations (the control parameters are the excitation emerge from the segment marked by P3B3 (period: 3, amplitudes). Next, the details are demonstrated through number of the branches: 3). Similarly, three branches a specific example with frequencies ω = 3 and appear from the segments P4 B3 and P3B3 on the hor- R1 ω = 2 (see again Fig. 1 for the driving signals). izontal axis, and two branches are merged together at R2 It must be emphasized that the amplitudes for both fre- segments P2B2 and P3B2 on the vertical axis in the quencies of the driving can be arbitrary; that is, there last row and second column of Fig. 2, better to be seen 123 Nonlinear oscillators by dual-frequency driving Fig. 2 List of high-resolution bi-parametric plots as a function of the pressure amplitudes at different relative-frequency pairs organized in rows. Left column: Lyapunov- exponent diagrams separating the chaotic and periodic solutions. Right column: Periodicity diagrams, where periodic solutions are plotted with different colours up to period-6. Period-7 or higher solutions (including chaos) can be found in the black regions. The pressure amplitudes P and P are A1 A2 given in bar in the enlargement and extension in Fig. 3. Observe In order to get a deeper insight into the branch- also that the number of the branches is exactly the ing mechanism at the relative-frequency pair ω = R1 same as the value of the corresponding relative fre- 3,ω = 2, computations were repeated with much R2 quency along the respective axis. This phenomenon higher resolution of the pressure amplitudes (Fig. 3). will be referred to as branching mechanism through- The ranges of the control parameters are P ∈ A1 out the paper and will be investigated in more detail for [0, 8] (bar) and P ∈[0, 5] (bar) with resolutions A2 the case ω = 3,ω = 2. 4001 and 2501, respectively. The number of initial R1 R2 123 F. Hegedus ˝ et al. Fig. 3 Bi-parametric periodicity diagram at the relative- grid in the y – y phase plane. The colour code is the same as in 1 2 frequency pair ω = 3,ω = 2 with increased resolution the right column of Fig. 2. The pressure amplitudes P and P R1 R2 A1 A2 (4001 × 2501) of the control parameters P and P . The num- are given in bar. (Color figure online) A1 A2 ber of the initial conditions is 5×5 = 25 defined on an equidistant conditions at a control parameter pair ( P , P ) is bordered each by a yellow (period-4) band at the right A1 A2 increased from 10 to 5 × 5 = 25 defined on an equidis- border (a period-doubled zone) followed by a thin black tant grid on the y – y phase plane. Thus, 2501×4001× line (indicating period-8 and higher periods including 1 2 25 ≈ 250 million initial-value problems have to be chaos that are encoded that way). Similarly, the three solved. This huge amount of computations is divided dark blue P3 bands are bordered each by a light blue into 5 × 8 = 40 (Δ P = 1 bar) segments and dis- (period-6) band, again indicating the start of a period- A1,2 tributed to different GPUs. doubling sequence. In addition, the three yellow P4 At first sight, Fig. 3 contains only slightly more bands are bordered by a thin black line, presumably information about the branching mechanism than the also the beginning of a period-doubling sequence to be bottom right diagram in Fig. 2, since it is hard to rep- seen only at higher resolution. Note that also other dark resent the many coexisting attractors in a single plot blue areas are bordered on one side by a period-doubled obtained by the increased number of initial conditions zone of light blue followed by black. Thus, the branches (in case of coexistence, the solution with the high- themselves are period-doubling objects. It should also est period is plotted). Investigating the system period- be noted that in any case, two of the three branches in by-period, however, this high-resolution bi-parametric each Pn B3, n = 2, 3, 4, turn to the left (lower driving scan becomes very helpful to understand the branching amplitudes P ) and one to the right (higher driving A1 mechanism (explored rigorously in the next sections); amplitudes P ). A1 in particular, to identify the connected branches with different periods confidently, such as the ones originat- ing from P4B3, P3B3, P2B3, P3B2 or P2B2, where, 5.2 Coexisting period-1 solutions for instance, P2B3 and P2B2 are connected. Typical colour combinations can be seen at the bor- Although the red period-1 domain in Fig. 3 shows no der of the branch families. The three green P2 bands are sign of multiple branches, the examination of its inner 123 Nonlinear oscillators by dual-frequency driving structure—the organization of the coexisting period- 1 orbits—reveals important features of the branching mechanism. They turn out to be general for other peri- ods as well. By filtering the results for period-1 orbits and by counting the number of the different attractors at a given control parameter pair, the structure of the coexisting period-1 solutions can be made visible in the P – P plane, see Fig. 4 with a magnification of A1 A2 the zone near the origin in the bottom panel. Domains with different colours represent areas with a differ- ent number of coexisting attractors up to four. These domains are bounded by saddle-node (SN, solid) and period-doubling (PD, dashed) curves computed by the boundary-value problem solver AUTO introduced in more detail later. The boundaries computed by AUTO are continuous curves. Some fractal-like shapes in the colour-coded domains are an artefact of the lack of a sufficient number of initial conditions or the presence of extremely long transients during the initial-value- problem computation Crossing one of the boundaries means the increment or decrement of the number of the attractors by one. Fig. 4 Number of coexisting period-1 orbits in the P – P A1 A2 Multiple stable period-1 orbits in a bi-parametric parameter plane at the relative-frequency pair ω = 3,ω = R1 R2 space usually originate via cusp bifurcations forming 2. The saddle-node (SN, solid) and period-doubling (PD, dashed) pairs of SN points (hysteresis) and giving rise to hys- curves, computed by the boundary-value problem solver AUTO, are the boundaries of the domains with a different number of teresis zones. The overlapping of such hysteresis zones coexisting attractors. The pressure amplitudes P and P are A1 A2 may produce domains with a high number of coexist- giveninbar ing period-1 attractors [53]. The sole cusp bifurcation in the bottom panel of Fig. 4, however, cannot explain the mosaic-type pattern of the period-1 structure. The cation points are indicated by PF, SN and PD, respec- investigation of the two special cases, where one of the tively. amplitudes of the dual-frequency driving P or P The period-reducing phenomenon discussed in detail A1 A2 is zero (vertical or horizontal axis in Fig. 4), can help in Sect. 4 plays an important role in the period-1 bifur- to understand this pattern. cation structure. Observe that the applied frequency combination is the same as in Fig. 1; furthermore, the first component of the dual-frequency driving is zero. Consequently, the driving is a purely harmonic 5.3 The limiting case P = 0 and the function proportional to the blue curve in Fig. 1, and A1 period-reducing phenomenon an originally period-2k solution is observed only as period-k orbit. Thus, the first bifurcation of the sta- The bifurcation structure with P = 0 bar is shown in ble period-1 curve originating from the equilibrium A1 Fig. 5, where the first component of the Poincaré sec- solution y (both amplitudes are zero) is a pitchfork 1 E tion Π( y ) is plotted as a function of P . The colour point PF, instead of a period-doubling point PD, which 1 A2 code is the same as in the cases of Fig. 3 and the right gives birth of two period-1 branches 2 × P1. Only the column of Fig. 2. For instance, red, green and blue subsequent bifurcations exhibit a Feigenbaum cascade curves are period-1, -2 and -3 orbits, respectively. Some (2 × P2 → 2 × P4 → ··· ). Similarly, the first bifur- of the periodic orbits are also marked by N ×PX, where cation point in a cascade must be a symmetry-breaking N is the number of the coexisting attractors of period X. bifurcation in case of a symmetric orbit in a symmetric The pitchfork, saddle-node and period-doubling bifur- nonlinear oscillator [54]. 123 F. Hegedus ˝ et al. Fig. 5 One-dimensional bifurcation structure with P = 0at periodic orbits are marked by N × PX,where N is the number of A1 the relative-frequency pair ω = 3, ω = 2, where the first the coexisting attractors of period X. The pitchfork, saddle-node R1 R2 component of the Poincaré section Π( y ) is plotted as a function and period-doubling bifurcation points are marked by PF, SN of the control parameter P . The colour code is the same as in and PD, respectively A2 the cases of Fig. 3 and the right column of Fig. 2. The relevant 1.4 Two examples of bubble radius versus time curves 1.3 y (τ ) are plotted in the top panel of Fig. 6 so as to 1.2 help understanding the period-reducing mechanism. The blue and black period-1 solutions are initiated 1.1 at P = 1 bar from the black dots marked by the A2 bold numbers 1 and 2 in Fig. 5, respectively. The blue 0.9 harmonic function in the bottom panel in Fig. 6 is 0.8 the purely harmonic driving of the system, see also 0.7 Fig. 1. It is clear that if the Poincaré sections are taken at time instances τ = n · T (n = 0, 1, 2,...), 0.6 0 0.5 1 1.5 2 2.5 3 which is used during the present simulations, the solu- [-] tions are regarded as period-1 orbits. In contrast, if the =2 R2 Poincaré sections are taken at time instances τ = n · T 1 n 2 (n = 0, 1, 2,...), which is a common choice in case of a purely harmonic driving, the solutions are period-2 -1 orbits. To prevent ambiguity in the definition of the peri- odicities, solutions will be referred to as “originally” -2 0 0.5 1 1.5 2 2.5 3 period- N orbits if the latter definition of the Poincaré [-] section is applied. Although the blue and black attrac- tors in Fig. 6 are different (each has its own basin of Fig. 6 The period-reducing phenomenon. Top panel: Bubble radius versus time curves y (τ ) (pressure amplitude P = 0bar 1 A1 attraction), the only difference between them is a phase P = 1 bar) starting at the black dots at τ = 0marked by A2 shift in time by one period of the driving (Δτ = 1.5). the bold numbers 1 and 2 in Fig. 5. Bottom panel: The sig- Altogether two types of bifurcation scenarios of nal of the single-frequency driving with period τ = T .The periodic orbits are observable in Fig. 5. For odd basic global Poincaré section applied during the present computations is defined by T (dashed, red vertical line). (Color figure online) periods (period before a bifurcation cascade takes place), the first point is a pitchfork bifurcation followed by a Feigenbaum period-doubling cascade (PF → y [-] p [bar] A2 Nonlinear oscillators by dual-frequency driving (A) (B) 2 0.7 =3 R1 0.6 =2 R2 1 1.5 P1 2 0.5 P =0 bar E 1 A1 PD P1 +P1 1 2 0.4 P =0.3 bar 0.3 E,u A1 E SN P =0.6 bar P1 P1 A1 0.2 0.5 PD P1 PF 0.1 E,u 2 P1 +P1 P =0 bar 1 2 A1 -0.1 2E -0.5 -0.2 012345 0.15 0.2 0.25 P [bar] P [bar] A2 A2 Fig. 7 Period-1 bifurcation curves computed by the boundary- node (SN) bifurcation points, respectively. The left panel (A) value problem solver AUTO. The black and red curves are stable shows the case with P = 0, while the right panel (B) A1 and unstable solutions, respectively. The asterisk, crosses and presents results at pressure amplitude P = 0.3 bar and at A1 dots are the pitchfork (PF), period-doubling (PD) and saddle- P = 0.6 bar. (Color figure online) A1 PD → PD → ··· ). Some examples are the already In order to illustrate this effect, it is useful to solve investigated period-1 curve P1 starting from y ,the and compute complete bifurcation curves of periodic 1 E period-3 orbits P3 originating from a SN point at orbits with a boundary-value-problem solver. As an P = 1.5 bar or the period-5 branches P5 appearing example, in Fig. 7A, the second component of the A2 again via a SN point at P = 2 bar. Observe that the Poincaré section Π( y ) of the period-1 curves is plot- A2 2 PF points have no influence to the colour, i.e. period. ted as a function of the control parameter P com- A2 In case of even basic periods, the period is immedi- puted by the bifurcation analysis and continuation soft- ately halved (period-reducing phenomenon), and the ware AUTO [36]. Here, the first component of the bifurcation scenario is a classic Feigenbaum cascade pressure amplitude P is still zero. AUTO is capable A1 without an initial PF point (PD → PD → ··· ). An of tracking down whole bifurcation curves including example here is the originally period-6 solution 2 × P3 the unstable part even if they contain multiple turn- generated at P = 2.5 bar. ing points, and it can detect the bifurcations and their A2 Through the discussion of Figs. 5 and 6, it has been types. This is the reason why AUTO is commonly used shown that the definition of the Poincaré section has an to study the bifurcation structure of nonlinear systems important influence on the apparent periodicity of the [5,53,55–61]. In Fig. 7, the stable and unstable parts of solutions. However, problems only appear in the two the period-1 branches are indicated by the black and limiting cases of either P = 0 bar or P = 0 bar. In red curves, respectively. The curve originating from A1 A2 the first special case of P = 0 bar, a suitable choice y becomes unstable via a pitchfork bifurcation PF A1 2 E can be either T or T , see the bottom panel of Fig. 6. at P = 0.19 bar producing two stable coexisting 2 A2 Since the addition of any small value to the first pres- period-1 branches. These stable curves lose their stabil- sure amplitude P will destroy the special structure ity approximately at P = 3 bar by a period-doubling A1 A2 introduced in Fig. 5, T will be no longer a period of point PD. the dual-frequency driving but only T . Therefore, in order to keep the definition unified, the usage of T as 5.4 Perturbation of the limiting case P = 0 A1 a global Poincaré section is advantageous. This, at first sight only technical, extension of the definition to the The breakup of the structurally unstable pitchfork bifur- limiting, purely harmonic-driving case is of great help cation due to the effect of a nonzero value of the ampli- in ordering and explaining the bifurcation scenarios of tude of the first harmonic component of the excitation dual-frequency driven systems. ( P = 0.3 bar) is clearly seen in Fig. 7B. The curves A1 E 1 marked by P1 and P1 in Fig. 7A form a single curve 1 2 (y ) [-] (y ) [-] 2 F. Hegedus ˝ et al. Fig. 8 One-dimensional bifurcation structure with P = 0at periodic orbits are marked by N × PX,where N is the number A2 the relative-frequency pair ω = 3, ω = 2, where the first of the coexisting attractors of period X. The saddle-node and R1 R2 component of the Poincaré section Π( y ) is plotted as a function period-doubling bifurcation points are marked by SN and PD, of the control parameter P . The colour code is the same as in respectively. (Color figure online) A1 the cases of Fig. 3 and the right column of Fig. 2. The relevant E ,u E 1 P1 + P1 in Fig. 7B. Likewise, the curves P1 and lower periodicities, and how these multiple orbits are 1 2 1 P1 presented in Fig. 7A form again a single curve separated under the perturbation of the amplitude of the E ,u other harmonic component. In order to understand the P1 + P1 in Fig. 7B composed by a stable and an 1 2 global picture, however, the investigation of the other unstable part separated by a saddle-node bifurcation special case, where the control parameter is P and (black dot) at the turning point. Elevating P from A1 A1 the amplitude of the second component P is zero, 0.3 bar to 0.6 bar, a pair of SN points appears along the A2 E 1 is necessary. The corresponding bifurcation diagram branch P1 + P1 forming a hysteresis. Observe that 1 2 E ,u E 1 2 is shown in Fig. 8, where the first component of the both curves P1 + P1 and P1 + P1 have their 1 2 1 2 Poincaré section Π( y ) is plotted as a function of the own separate evolution as the amplitude P changes; A1 pressure amplitude P . The colour code and the nota- A1 consequently, the black and blue curves in the top of tion system are the same as in Fig. 5.Atfirstsight—in Fig. 6 must become different (not only by a phase shift view of (2) with θ = 0—both cases ( P = 0 and A1 in time) for P > 0 bar as they have their own sepa- A1 P = 0) seem exchangeable. However, switching to A1 rate evolution as well. All the pitchfork structures in the the other frequency component, different bifurcation bifurcation scenarios initiated from odd basic periods scenarios develop (if the frequencies are distinct). are broken up in the same way as presented in Fig. 7B, There are two types of bifurcation scenarios of peri- except for the formation of the hysteresis which not odic orbits. If the basic (original) period is divisible necessarily happens. by 3, then the period reducing takes place immediately decreasing the period of the orbits from 3k to k followed by a Feigenbaum cascade as usual (PD → PD → ··· ). 5.5 The limiting case P = 0 A2 Notable examples are the branches 3 × P1 and 3 × P3 both appearing via SN bifurcations at P = 1.2 bar A1 Shortly in this section, the complex development of and at P = 4.1 bar, respectively. For other basic A1 the period-1 structure in a wide range of the amplitudes periods, there is no period-reducing phenomenon, but P and P will be discussed in detail. Up to now, it is A1 A2 only a period-doubling cascade. Examples of this sce- important only to understand how the period-reducing nario are the period-1 curve departing from the equi- mechanism works for single-frequency driving, how it librium solution y , the period-4 orbits generated at 1 E decomposes periodic solutions into multiple orbits of 123 Nonlinear oscillators by dual-frequency driving 1.4 stable up to the pressure amplitude P ≈ 7 bar. It A1 1.3 loses its stability there via a period-doubling bifur- cation. The originally period-3 family of solutions 1.2 3 × P1 appearing through a SN point is decoupled into 1.1 three period-1 branches existing between the ampli- tudes P ≈ 1.19 bar and P ≈ 4.67 bar. Therefore, A1 A1 0.9 there is a wide range of amplitudes P where four A1 period-1 attractors exist. Figure 9 shows these coexist- 0.8 ing orbits at P = 2 bar starting from the black dots A1 0.7 in Fig. 8. The dashed curve is the period-1 orbit related 0.6 to the curve originating from y . The black, red and 1 E 0 0.5 1 1.5 2 2.5 3 [-] blue solid curves are the orbits of the period-reduced 3 × P1 branches, where the only difference between Fig. 9 The period-reducing phenomenon. Dimensionless bub- them is again a phase shift in time, in this case of one ble radius versus time curves of four coexisting stable period-1 or two periods T (Δτ = 1or2). orbits at pressure amplitude P = 2 bar starting from the black A1 1 dots shown in Fig. 8 at τ = 0. Between the black, red and blue solid curves there is only a phase shift in time with Δτ = 1 or 2. The dashed curve is a degenerate (i.e. repeating only after T instead of after T ), small-bubble-radius-amplitude case that 1 5.6 Perturbation of the limiting case P = 0 A2 exists stably up to P ≈ 7 bar where the oscillation period A1 doubles. (Color figure online) Again, the boundary-value-problem solver AUTO can help to understand the behaviour of the system under a small perturbation of the amplitude of the second, P = 1.5 bar or the period-5 branches initiated at harmonic component of the driving ( P > 0). The A1 A2 P = 4.6 bar. Observe that there are no pitchfork or bifurcation curves related to the four period-1 orbits A1 other special bifurcations present in these cases. can be seen in Fig. 10Afor P = 0 bar. The black A2 Figure 8 explains how the coexistence is possible of and red curves mean stable and unstable orbits, respec- altogether four period-1 stable solutions highlighted in tively. The bifurcation of the saddle-node (SN) and Fig. 4 by the yellow domain. The period-1 orbit P1 period-doubling (PD) points are also detected. They emerges from the equilibrium point y and remains are marked by dots (SN) and crosses (PD). Observe 1 E 2.5 1.4 (B) =3 (A) P =0.1 bar R1 A2 1.2 PD =2 2 1 R2 P1 3 1 P =0 bar A2 1.5 0.8 0.6 SN 0.4 P1 0.5 0.2 P1 3 0 2 -0.2 P1 2E -0.5 -0.4 01234567 8 0.5 1 1.5 2 P [bar] P [bar] A1 A1 Fig. 10 Period-1 bifurcation curves computed by the boundary- while the right panel (B) presents results of the 3 × P1curves value problem solver AUTO. The black and red curves are stable at the pressure amplitude of P = 0.1 bar. The arrows in panel A2 and unstable solutions, respectively. The crosses and dots are the (B) show the motion of the SN bifurcation points, two of them period-doubling (PD) and saddle-node (SN) bifurcation points, towards lower P , one of them towards larger P . (Color figure A1 A1 respectively. The left panel (A) shows the case with P = 0, online) A2 (y ) [-] y [-] 2 1 (y ) [-] 2 F. Hegedus ˝ et al. that the P -ranges of the stable (black) segments in A1 Fig. 10A coincide with the domain of existence of the corresponding parts of the red curves in Fig. 8.Onthe right hand side of Fig. 10, only the evolution of the curves 3 × P1 is presented by increasing the ampli- tude P fromzeroto0.1 bar. It is clear that while the A2 saddle-node bifurcations (black dots) take place at the same value of P when P = 0 bar (see the vertical A1 A2 thin line in Fig. 10B), the locations of these bifurcation points differ when P > 0. This means that each of A2 the three bifurcation curves have their own evolution as the amplitude P changes. Consequently, the black, A2 red and blue solid curves in Fig. 9 become different (not Fig. 11 3D representation of the stable period-1 orbits (Poincaré only by a phase shift in time) for P > 0. A2 section points of the bubble wall velocity) above the parameter plane P – P by filtering the results presented in Fig. 3 by A1 A2 period, here period-1. The period-1 orbits can be separated into four different layers marked by the blue, black, red and green surfaces. The indication of the curves in the Π( y )– P (Fig. 10) 2 A1 6 Period 1: the complete picture m or Π( y )– P (Fig. 7)planesis P1 ,where n is the original 2 A2 period and m = 1, ..., n is a serial number, P1 means period- 1 solution. The black surface originating from the equilibrium Due to the exhaustive discussion in the preceding sec- E 1 solution of the system is marked by P1 (instead of P1 to 1 1 tions, it is clear now how the system behaves under show the connection with the equilibrium point). (Color figure single-frequency driving when a second, harmonic online) component is added as a perturbation. The huge amount of results presented in Fig. 3 can help to establish a global overview of the evolution of all the period-1 Figure 11 indicates a very interesting phenomenon, solutions by means of filtering the results by period, namely the originally period-2 solutions on P1 can here period-1. Figure 11 shows these period-1 solu- be transformed into the originally period-3 solutions tions represented in a three-dimensional plot, where the on P1 through the blue surface. Similarly, orbits on 2 2 second component of the Poincaré section, Π( y ),the P1 can be transformed into orbits on P1 via the 2 3 bubble wall velocity, is plotted as a function of the two green surface. These are smooth transformations by main control parameters P and P . The different means of the continuous change of the amplitudes P A1 A2 A1 period-1 orbits can be decomposed into several layers and P of the dual-frequency driving. Naturally, all A2 indicated by the blue, black, red and green surfaces. these orbits are considered as period-1 attractors due In the limit cases when one of the pressure amplitudes to the specific choice of the global Poincaré section. is zero (planes Π( y )– P or Π( y )– P ), the points Consider, however, the following scenario of parame- 2 A1 2 A2 are plotted with bigger markers to be able to distin- ter variations visualized by the supplementary movie guish such special solutions more easily. These solu- file (Online Resource 1) where the top panel shows the tions, which form complete bifurcations curves pre- dimensionless bubble radius curves y as a function of sented in Figs. 7A and 10A, are marked by P1 , where the dimensionless time τ , and the bottom panel repre- n is the original period and m = 1,..., n is a serial sents the dual-frequency driving signal P (τ ) (exclud- number, P1 means period-1 solution. When both P ing the static ambient pressure). A1 and P are zero then the system is in equilibrium and Firstly, let us start the investigation with a solu- A2 the dimensionless bubble wall velocity and its Poincaré tion lying on the bifurcation curve 2 × P1inFig. 5 1 2 section y are zero (it is the origin of the diagram). ( P1 or P1 in Fig. 11) so that the amplitude of 2 E 2 2 The black period-1 surface originating from this point is the single-frequency excitation is somewhere between E 1 marked by P1 (instead of P1 , emphasizing the equi- P = 1.4 bar and P = 2.9 bar ( P = 0 bar) with A2 A2 A1 1 1 librium origin E). There is only one curve of this type, relative frequency ω = 2. Specifically, in Online R2 which is a surface of small-amplitude bubble oscilla- Resource 1, the pressure amplitude is chosen to be tions. P = 2.5 bar. In this video, besides the originally A2 123 Nonlinear oscillators by dual-frequency driving Fig. 12 A 3D representation of the bifurcation structure, com- limit cases of P = 0bar or P = 0 bar. Panel B is an enlarge- A1 A2 puted by AUTO, of the blue and black layers presented in Fig. 11. ment of the near-origin zone to better show the bifurcation-point The black and red curves are the stable and unstable solutions, area. (Color figure online) respectively. The blue and green curves form the skeleton of the period-2 (red) attractor, which is the initial state of the tinuous domain in the pressure amplitude—frequency transformation, an originally period-3 (blue) attractor parameter plane [53,61]. Therefore, a smooth transfor- coexists regarded as a final target of the control prob- mation of a period-2 solution into a period-3 orbit of a lem. The colour code of these attractors is the same single-frequency driven system is possible by a tempo- both in Fig. 5 and in the top panel of Online Resource rary addition of a second harmonic component to the 1; in addition, by comparing the top and bottom panels, driving. This is an efficient way to control multistability the periodicities of the attractors are clearly visible. in harmonically driven nonlinear oscillators. Accord- Secondly, guide the period-2 orbit smoothly to the ing to the best of the authors’ knowledge, this phe- 1 2 curve 3 × P1 shown in Fig. 8 (either P1 or P1 in nomenon has not yet been published in the literature. 3 3 Fig. 11) by varying the pressure amplitudes continu- During this third step in Online Resource 1, the change ously. By the end of the second operation, the sys- of P is not necessary and kept constant. Observe how A1 tem is again driven by a single frequency (ω = 3) the initial and the final parameter set become identi- R1 with pressure amplitude between P = 1.2 bar and cal; that is, single-frequency driving with amplitude A1 P = 4.7 bar ( P = 0 bar). In Online Resource 1, P = 2.5 bar at ω = 2. The only difference is A1 A2 A1,2 R1,2 this transformation takes place along a circle in the the state of the system: operation on the period-3 (blue) P – P plane with radius 2.5 bar; that is, the final instead on the period-2 (red) attractor. A1 A2 value of the first pressure amplitude is P = 2.5 bar. In order for a better visualization of the surfaces pre- A1 The continuously changing trajectory during the trans- sented in Fig. 11, hundreds of bifurcation curves are formation is marked by the black curve in the top panel. computed by AUTO departing from the results shown Thirdly, keep the single-frequency driving but change in Fig. 10A. These surfaces, composed by lines, are both the amplitude P and the relative frequency ω summarized in Figs. 12, 13 and 14. The stable and A1 R1 of the driving (in general) to transform the solution unstable parts are the black and red curves, respectively. from the red curve 3 × P1inFig. 8 to the blue curve As a skeleton, all the bifurcation curves of the limit P3inFig. 5. This third step is possible, since both the cases ( P = 0 bar or P = 0 bar) are presented by A1 A2 3 × P1 and the P3 curves are related to the same sub- the blue (stable) and the green (unstable) curves, which harmonic resonance of order 1/3, which forms a con- can help to identify the surfaces more easily (compare 123 F. Hegedus ˝ et al. (A) (B) 2.5 1.5 1.5 1 0.5 0.5 -0.5 0 3 SN P1 5 PF 4 P1 8 8 1 5 1 4 0 0 P [bar] 0 P [bar] 1 A2 A2 P [bar] P [bar] A1 A1 Fig. 13 A 3D representation of the bifurcation structure, com- The blue and green curves form the skeleton of the limit cases puted by AUTO, of the red layer presented in Fig. 11. The black of P = 0bar or P = 0 bar. (Color figure online) A1 A2 and red curves are the stable and unstable solutions, respectively. (A) (B) 2.5 1.5 1.5 0.5 0.5 -0.5 P1 2 2 P1 4 -0.5 P1 P1 3 3 3 PF 3 SN 2 P1 0 0 0 P [bar] 0 P [bar] A2 P [bar] P [bar] A1 A1 A2 Fig. 14 A 3D representation of the bifurcation structure, com- The blue and green curves form the skeleton of the limit cases puted by AUTO, of the green layer presented in Fig. 11. The black of P = 0bar or P = 0 bar. (Color figure online) A1 A2 and red curves are the stable and unstable solutions, respectively. also with Figs. 7A, 10A). As usual, the pitchfork (PF), the right hand side is a magnification of the interesting saddle-node (SN) and period-doubling (PD) bifurca- parts. tion points are marked by asterisk, dots and crosses, From Fig. 12, it is clear that the blue and black sur- respectively. The left hand side of these figures repre- faces in Fig. 11 are connected for very small values sents the surfaces in the full parameter domain, while of P . The separation takes place approximately at A1 P = 0.33 bar via a cusp bifurcation forming a hys- A1 (y ) [-] (y ) [-] 2 2 (y ) [-] (y ) [-] 2 Nonlinear oscillators by dual-frequency driving panel of Fig. 15, there are two branches along which the following conversion occurs: 2 × P3 = P6 ↔ 3 × P3 = P9; that is, an originally period-6 solution is transformed into an originally period-9 orbit. Simi- larly, the bottom panel of Fig. 15 shows four bands con- verting originally period-8 orbits into period-12 orbits (2 × P4 = P8 ↔ 3 × P4 = P12). As a generaliza- tion, the present study found that applying the relative- frequency combination ω = 3 and ω = 2, a trans- R1 R2 formation between period-2k and period-3k attractors is possible, here k = 1, 2,.... 8 Winding numbers and compatibility of periodic orbits Naturally, not every kind of period-2k and period-3k orbits can be transformed into each other; in other words, there is a “compatibility” condition. Observe that the first PD points of the period-1 branches in Fig. 5 (2 × P1) and in Fig. 8 (3 × P1) are connected via the period-doubling curve PD shown in the top panel of Fig. 4. It is well known in the literature that the quan- tity called winding number w = nk/mk is invariant and Fig. 15 Number of the coexisting period-3 orbits (top panel) and period-4 orbits (bottom panel) in the P – P parameter plane A1 A2 rational near a period-doubling or a saddle-node point at relative frequencies ω = 3and ω = 2 R1 R2 [62,63]. Here, mk is the period of a given attractor and nk is the number of the rotations of a nearby trajectory around the attractor during one period. Therefore, to teresis, see also the bottom panel of Fig. 4. An example each bi-parametric PD or SN bifurcation curve a sin- for this hysteresis, which exists between the regions of gle, rational winding number can be associated. When a P > 0.33 bar and P < 1.19 bar, is already intro- A1 A1 transformation between period-2k and period-3k orbits duced in Fig. 7Bfor P = 0.6 bar. As expected, the A1 is possible, then it must be possible between their cor- red, standalone surface in Fig. 11 has no connections responding bifurcation points as well, see, for instance, with other surfaces, see Fig. 13. The green surface in the PD curve in Fig. 4 which has the winding number Fig. 11 together with its unstable counterpart forms a w = 1/2. Consequently, the transformation is possible simple banded, dual-layer surface presented in Fig. 14. only if the winding numbers of such bifurcation points are compatible; that is, when they are equal. 7 Higher periods 9 Discussion If the smooth transformation works between certain period-2 and period-3 orbits, there must be other com- According to the literature, the only option to directly binations of periodic solutions where the transforma- target an attractor in case of multistability is to apply a tion is possible. In order to investigate these other pos- feedback control technique. The main drawback of the sibilities, the number of the coexisting period-3 and feedback control methods is that explicit knowledge of period-4 attractors is plotted in the top and the bottom the state space (modified targeting method [33], bush- panels of Fig. 15, respectively. The colour codes are the like paths to a pre-selected attractor [64], reinforcement same as in case of Fig. 4 (period-1 orbits). Let us focus learning [65]) or the properties of the targeted attrac- only on those parts of the complex bifurcation struc- tor (trajectory selection by a periodic feedback [34]) is tures where the transformations take place. In the top necessary. Control of multistability is possible by delay 123 F. Hegedus ˝ et al. feedback as well, which can be combined with chaos cation structure of many nonlinear oscillators have the control easily [66]. In that case, the feedback stabi- same basic building blocks [70]. If there is no informa- lizes one of the unstable orbits embedded in a chaotic tion about a state variable, the orientation in the topol- attractor. Unfortunately, it just stabilizes an unstable ogy can be supported; for instance, by measuring and orbit and does not target directly an already existing Fourier transforming an indirect quantity connected to attractor. Parenthetically, the feedback control can sta- the dynamics (in case of bubbles their emitted pres- bilize an attractor against external noise by employing sure signal). The peaks presented in the spectra can the Jacobian of the attractor as a feedback [67]. How- provide valuable information about the nature of the ever, the state of the system must be already near to actual attractor [71], see also the studies of Lauterborn the desired orbit; that is, the application of one of the and Cramer [72] and Lauterborn and Suchla [73] where above targeting techniques is mandatory. complete period-doubling scenarios could be measured If feedback control is possible then the aforemen- in this way of harmonically forced bubbles represented tioned controls work very well; otherwise, the appli- via spectral bifurcation diagrams. Another requirement cation of a non-feedback method (including a stochas- related to the topology is that the transformed trajecto- tic one) is inevitable. The main advantage of the non- ries must be compatible in terms of winding numbers feedback control technique presented in this study is as already discussed in Sect. 8. the ability to target the desired attractor directly, which It must be emphasized that it is possible to combine is not possible with traditional non-feedback controls. the present non-feedback control technique with other For instance, the attractor selection by short pulses [31] methods. For instance, with stochastic control, attrac- (kick the system) introduces randomness to the con- tors with small basins can be annihilated to “simplify” trol; that is, it cannot be guaranteed that the final state the state spaces during the transformation. Moreover, of the system after the kick is the desired attractor. even if feedback control is possible in an application, Moreover, extremely long transients in the presence the method introduced here is still a good candidate for of transient chaos [68] are also possible. With the con- targeting a desired attractor, and after the continuous tinuous transformation of attractors such transients can transformation, the attractor can be stabilized against be avoided. The attractor annihilation techniques by external noise with the feedback control. Similarly, pseudo-periodic forcing [69] or by harmonic pertur- the delay feedback control can stabilize an unstable bation of a state variably or a parameter [32] cannot orbit that replaces the originally stable chaotic attrac- target an attractor directly either. Roughly speaking, tor. After that, the continuous transformation technique both techniques destruct attractors in the order of the can be applied between the resulting coexisting stable size of their basin of attraction with increasing forcing states if necessary and if it is possible. In this way, chaos or perturbation amplitudes. control and targeting an attractor directly is accom- The applicability of the method presented here has plished simultaneously. There are case studies where also requirements/limitations. Firstly, in many cases, the application of two different kinds of controls may there is the requirement that the system should not be be successful when the usage of only one of the control modified appreciably during the control. Therefore, all fails [74]. non-feedback techniques apply only small perturba- Since the control technique is demonstrated numer- tions. It is absolutely not fulfilled in the present method ically only on a specific system (Keller–Miksis equa- as large deviations from the original system parameters tion describing bubble dynamics) and on a very spe- is usually necessary. However, by the end of the control, cific example, many questions arise: Is it possible to the system goes back exactly to the same parameter set “bypass” somehow the winding number compatibility and additional control is no longer required to sustain condition? Is it possible to transform a chaotic attrac- the desired state (attractor). It must be also noted that if tor into a stable periodic orbit and vice versa? Is this the basin of attraction of the targeted attractor is small method applicable to other systems? In order to answer then slow transformation may be necessary. these question correctly, additional detailed analyses Secondly, due to the large parameter variation during are mandatory; however, preliminary conjectures can the control, detailed knowledge on the topology of the be made. periodic orbits in a large parameter space must be avail- Even if a continuous transformation between two able. This is usually a minor problem since the bifur- coexisting attractors does not exist applying a single- 123 Nonlinear oscillators by dual-frequency driving frequency pair, it may be possible to find a transfor- that two attractors which cannot be exchanged by tun- mation by using subsequent frequency pairs one after ing the system parameters may be transformed contin- another, and reaching the desired attractor step-by-step uously into each other by temporarily adding a second hopping onto “intermediate” attractors. The transfor- sinusoidal component to the driving. The process to mation of chaotic attractors into periodic orbits fits into target a desired attractor directly was tested at a rele- the large topic of chaos control [21]. Since chaos con- vant example, the Keller–Miksis equation—a second- trol is possible with the addition of a second sinusoidal order ordinary differential equation describing bubble component to the driving [24,75], the main issue is that dynamics—and visualized via an animation (Online how chaos is annihilated (a smooth inverse Feigenbaum Resource 1). If this technique is proven to be applicable cascade or a sudden crises) during the parameter tuning, to other systems as well, it may become an excellent and what is the resulting periodic orbit which is inter- candidate to accomplish non-feedback control without esting from the winding number compatibility point of long transients and/or randomness. view. It is well known that chaos control is possible Acknowledgements Open access funding provided by Max merely by adjusting the phase shift between the sinu- Planck Society. This paper was supported by the János Bolyai soidal components of the driving [76]; therefore, the Research Scholarship of the Hungarian Academy of Sciences effect of the phase shift must also be thoroughly inves- and the Deutsche Forschungsgemeinschaft (DFG) Grant No. ME tigated. It can also play an important role between the 1645/7-1. transformation of the period-2 and period-3 orbits itself Compliance with ethical standards shown in Fig. 11. The period-2 solutions can always be transformed into one of the period-3 orbits through the Conflict of interest The authors declare that they have no con- blue or the green surfaces. In contrast, if the period-3 flict of interest. attractor is in a phase in time described by the red sur- Open Access This article is distributed under the terms of the face then transformation is not possible in the opposite Creative Commons Attribution 4.0 International License (http:// direction (from period-3 to period-2), see again Fig. 11. creativecommons.org/licenses/by/4.0/), which permits unrestricted It is very likely that the control technique of con- use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, tinuous transformation is possible for other systems provide a link to the Creative Commons license, and indicate if as well. As it was already mentioned, the topology of changes were made. many nonlinear oscillators is composed by the same building blocks. Moreover, in recent publications, a phenomenon called replication of periodic domains is References reported by applying a second sinusoidal component to the driving [40,77], or a single periodic perturbation in 1. Pisarchik, A.N., Feudel, U.: Control of multistability. Phys. case of maps [78,79]. We believe that the replication of Rep. 540(4), 167 (2014) 2. Angeli, D., Ferrell, J.E., Sontag, E.D.: Detection of multista- periodic domains and the branching mechanism intro- bility, bifurcations, and hysteresis in a large class of biolog- duced in Sect. 5.1 have the same dynamical origin. ical positive-feedback systems. Proc. Natl. Acad. Sci. USA Finally, it is worth mentioning that any application 101(7), 1822 (2004) using dual-frequency driving can benefit (at least indi- 3. Ullner, E., Koseska, A., Jürgen, K., Volkov, E., Kantz, H., García-Ojalvo, J.: Multistability of synthetic genetic net- rectly) from the results presented in this study due to works with repressive cell-to-cell communication. Phys. the very detailed investigation. Some possible appli- Rev. E 78, 031904 (2008) cations are acoustic cavitation [80–83], Faraday waves 4. Shiau, Y.H., Peng, Y.F., Hwang, R.R., Hu, C.K.: Multista- [84], stability of traveling beams [85,86]orlaser-driven bility and symmetry breaking in the two-dimensional flow around a square cylinder. Phys. Rev. E 60, 6188 (1999) dissociation of molecules [87]. 5. Hos, ˝ C.J., Champneys, A.R., Paul, K., McNeely, M.: Dynamic behaviour of direct spring loaded pressure relief valves in gas service: II reduced order modelling. J. Loss 10 Conclusion Prevent. Proc. 36, 1 (2015) 6. Ganapathisubramanian, N., Showalter, K.: Bistability, mushrooms, and isolas. J. Chem. Phys. 80(9), 4177 (1984) In the present study, a new non-feedback technique to 7. Yang, L., Dolnik, M., Zhabotinsky, A.M., Epstein, I.R.: control multistability in harmonically driven nonlinear Turing patterns beyond hexagons and stripes. Chaos 16(3), oscillators is introduced. It is based on the observation 037114 (2006) 123 F. Hegedus ˝ et al. 8. Braun, J., Mattia, M.: Attractors and noise: twin drivers of 30. Huisman, J., Weissing, F.: Biodiversity of plankton by decisions and multistability. Neuroimage 52(3), 740 (2010) species oscillations and chaos. Nature 402, 407–410 (1999) 9. Bathiany, S., Claussen, M., Fraedrich, K.: Implications of 31. Chizhevsky, V.N., Grigorieva, E.V., Kashchenko, S.A.: Opti- climate variability for the detection of multiple equilibria mal timing for targeting periodic orbits in a loss-driven CO2 and for rapid transitions in the atmosphere-vegetation sys- laser. Opt. Commun. 133(1), 189 (1997) tem. Clim. Dyn. 38(9), 1775 (2012) 32. Pisarchik, A.N., Goswami, B.K.: Annihilation of one of the 10. Sneppen, K., Mitarai, N.: Multistability with a metastable coexisting attractors in a bistable system. Phys. Rev. Lett. mixed state. Phys. Rev. Lett. 109, 100602 (2012) 84, 1423 (2000) 11. Lauterborn, W., Kurz, T.: Physics of bubble oscillations. 33. Shinbrot, T., Ott, E., Grebogi, C., Yorke, J.A.: Using chaos Rep. Prog. Phys. 73(10), 106501 (2010) to direct trajectories to targets. Phys. Rev. Lett. 65, 3215 12. Lieberman, M.A., Tsang, K.Y.: Transient chaos in dissipa- (1990) tively perturbed, near-integrable hamiltonian systems. Phys. 34. Jiang, Y.: Trajectory selection in multistable systems using Rev. Lett. 55, 908 (1985) periodic drivings. Phys. Lett. A 264(1), 22 (1999) 13. Kolmogorov, A.N.: On the persistence of conditionally peri- 35. Kaneko, K.: Dominance of milnor attractors and noise- odic motions under a small change of the Hamiltonian func- induced selection in a multiattractor system. Phys. Rev. Lett. tion. Dokl. Akad. Nauk. S.S.S.R. 98, 527 (1954). (Russian) 78, 2736 (1997) 14. Feudel, U., Grebogi, C., Hunt, B.R., Yorke, J.A.: Map with 36. Doedel, E.J., Oldeman, B.E., Champneys, A.R., Dercole, more than 100 coexisting low-period periodic attractors. F., Fairgrieve, T.F., Kuznetsov, Y.A., Paffenroth, R., Sand- Phys. Rev. E 54, 71 (1996) stede, B., Wang, X., Zhang, C.: AUTO-07P: Continuation 15. Gavrilov, N.K., Silnikov, L.P.: On three dimensional dynam- and Bifurcation Software for Ordinary Differential Equa- ical systems close to systems with structurally unstable tions. Concordia University, Montreal (2012) homoclinic curve. I. Math. U.S.S.R. Sbornik 17, 467 (1972) 37. Brennen, C.E.: Cavitation and Bubble Dynamics. Oxford 16. Newhouse, S.E.: Diffeomorphisms with infinitely many University Press, New York (1995) sinks. Topology 13(1), 9 (1974) 38. Niemeyer, K.E., Sung, C.J.: Accelerating moderately stiff 17. Wiesenfeld, K., Hadley, P.: Attractor crowding in oscillator chemical kinetics in reactive-flow simulations using GPUs. arrays. Phys. Rev. Lett. 62, 1335 (1989) J. Comput. Phys. 256, 854 (2014) 18. Ikeda, K.: Multiple-valued stationary state and its instability 39. Bonatto, C., Gallas, J.A.C., Ueda, Y.: Chaotic phase similar- of the transmitted light by a ring cavity system. Opt. Com- ities and recurrences in a damped-driven Duffing oscillator. mun. 30(2), 257 (1979) Phys. Rev. E 77(2), 026217 (2008) 19. Pisarchik, A.N.: Dynamical tracking of unstable periodic 40. Medeiros, E.S., de Souza, S.L.T., Medrano-T, R.O., Caldas, orbits. Phys. Lett. A 242(3), 152 (1998) I.L.: Replicate periodic windows in the parameter space of 20. Campos-Mejía, A., Pisarchik, A.N., Arroyo-Almanza, D.A.: driven oscillators. Chaos Solitons Fract. 44(11), 982 (2011) Noise-induced on-off intermittency in mutually coupled 41. de Souza, S.L.T., Lima, A.A., Caldas, I.L., Medrano-T, R.O., semiconductor lasers. Chaos Soliton. Fract. 54, 96 (2013) Guimaraes-Filho, ˝ Z.O.: Self-similarities of periodic struc- 21. Schöll, E., Schuster, H.G.: Handbook of Chaos Control. tures for a discrete model of a two-gene system. Phys. Lett. Wiley, New York (2008) A 376(15), 1290 (2012) 22. Mettin, R., Kurz, T.: Optimized periodic control of chaotic 42. Francke, R.E., Pöschel, T., Gallas, J.A.C.: Zig–zag networks systems. Phys. Lett. A 206(5), 331 (1995) of self-excited periodic oscillations in a tunnel diode and a 23. Mettin, R., Hübler, A., Scheeline, A., Lauterborn, W.: Para- fiber-ring laser. Phys. Rev. E 87(4), 042907 (2013) metric entrainment control of chaotic systems. Phys. Rev. E 43. Medrano-T, R.O., Rocha, R.: The negative side of Chua’s 51, 4065 (1995) circuit parameter space: stability analysis, period-adding, 24. Behnia, S., Sojahrood, A.J., Soltanpoor, W., Jahanbakhsh, basin of attraction metamorphoses, and experimental inves- O.: Suppressing chaotic oscillations of a spherical cavitation tigation. Int. J. Bifurcat. Chaos 24(09), 1430025 (2014) bubble through applying a periodic perturbation. Ultrason. 44. Celestino, A., Manchein, C., Albuquerque, H.A., Beims, Sonochem. 16(4), 502 (2009) M.W.: Stable structures in parameter space and optimal 25. Martínez, P.J., Euzzor, S., Gallas, J.A.C., Meucci, R., ratchet transport. Commun. Nonlinear Sci. Numer. Simul. Chacón, R.: Impulse-Induced Optimum Control of Chaos 19(1), 139 (2014) in Dissipative Driven Systems. ArXiv e-prints (2017) 45. Rech, P.C.: Period-adding structures in the parameter-space 26. Pyragas, K., Lange, F., Letz, T., Parisi, J., Kittel, A.: Stabi- of a driven Josephson junction. Int. J. Bifurc. Chaos 25(12), lization of an unstable steady state in intracavity frequency- 1530035 (2015) doubled lasers. Phys. Rev. E 61, 3721 (2000) 46. da Costa, D.R., Hansen, M., Guarise, G., Medrano-T, R.O., 27. Wieczorek, S., Krauskopf, B., Lenstra, D.: Mechanisms for Leonel, E.D.: The role of extreme orbits in the global organi- multistability in a semiconductor laser with optical injection. zation of periodic regions in parameter space for one dimen- Opt. Commun. 183(1), 215 (2000) sional maps. Phys. Lett. A 380(18), 1610 (2016) 28. Guevara, M., Glass, L., Shrier, A.: Phase locking, period- 47. Freire, J.G., Gallas, M.R., Gallas, J.A.C.: Stability mosaics doubling bifurcations, and irregular dynamics in period- in a forced Brusselator. Eur. Phys. J. Spec. Top. 226(9), 1987 ically stimulated cardiac cells. Science 214(4527), 1350 (2017) (1981) 48. Horstmann, A.C.C., Albuquerque, H.A., Manchein, C.: The 29. Li, E.: Chromatin modification and epigenetic reprogram- effect of temperature on generic stable periodic structures in ming in mammalian development. Nat. Rev. Genet. 3, 662– the parameter space of dissipative relativistic standard map. 673 (2002) Eur. Phys. J. B 90(5), 96 (2017) 123 Nonlinear oscillators by dual-frequency driving 49. Prants, F.G., Rech, P.C.: Complex dynamics of a three- 69. Pecora, L.M., Carroll, T.L.: Pseudoperiodic driving: elim- dimensional continuous-time autonomous system. Math. inating multiple domains of attraction using chaos. Phys. Comput. Simul. 136, 132 (2017) Rev. Lett. 67, 945 (1991) 50. Sewerin, F., Rigopoulos, S.: A methodology for the integra- 70. Scheffczyk, C., Parlitz, U., Kurz, T., Knop, W., Lauterborn, tion of stiff chemical kinetics on GPUs. Combust. Flame W.: Comparison of bifurcation structures of driven dissipa- 162(4), 1375 (2015) tive nonlinear oscillators. Phys. Rev. A 43(12), 6495 (1991) 51. Curtis, N.J., Niemeyer, K.E., Sung, C.J.: An Investigation of 71. Zhang, Y., Zhang, Y., Li, S.: Combination and simultaneous GPU-Based Stiff Chemical Kinetics Integration Methods, resonances of gas bubbles oscillating in liquids under dual- ArXiv e-prints (2016) frequency acoustic excitation. Ultrason. Sonochem. 35, 431 52. Stone, C.P., Alferman, A.T., Niemeyer, K.E.: Accelerating (2017) Finite-Rate Chemical Kinetics with Coprocessors: Com- 72. Lauterborn, W., Cramer, E.: Subharmonic route to chaos paring Vectorization Methods on GPUs, MICs, and CPUs, observed in acoustics. Phys. Rev. Lett. 47(20), 1445 (1981) ArXiv e-prints (2016) 73. Lauterborn, W., Suchla, E.: Bifurcation superstructure in a 53. Hegedus, ˝ F., Hos, ˝ C., Kullmann, L.: Stable period 1, 2 and model of acoustic turbulence. Phys. Rev. Lett. 53(24), 2304 3 structures of the harmonically excited Rayleigh–Plesset (1984) equation applying low ambient pressure. IMA J. Appl. Math. 74. Pisarchik, A.N., Martínez-Zérega, B.E.: Noise-induced 78(6), 1179 (2013) attractor annihilation in the delayed feedback logistic map. 54. Englisch, V., Parlitz, U., Lauterborn, W.: Comparison of Phys. Lett. A 377(42), 3016 (2013) winding-number sequences for symmetric and asymmetric 75. Zhang, Y., Zhang, Y.: Chaotic oscillations of gas bub- oscillatory systems. Phys. Rev. E 92(2), 022907 (2015) bles under dual-frequency acoustic excitation. Ultrason. 55. Fyrillas, M.M., Szeri, A.J.: Dissolution or growth of soluble Sonochem. 40(Part B), 151 (2018) spherical oscillating bubbles. J. Fluid Mech. 277, 381 (1994) 76. Yang, J., Qu, Z., Hu, G.: Duffing equation with two periodic 56. Hos, ˝ C., Champneys, A.R., Kullmann, L.: Bifurcation anal- forcings: the phase effect. Phys. Rev. E 53, 4402 (1996) ysis of surge and rotating stall in the Moore–Greitzer com- 77. Medeiros, E.S., de Souza, S.L.T., Medrano-T, R.O., Caldas, pression system. IMA J. Appl. Math. 68(2), 205 (2003) I.L.: Periodic window arising in the parameter space of an 57. Hegedus, ˝ F., Kullmann, L.: Basins of attraction in a har- impact oscillator. Phys. Lett. A 374(26), 2628 (2010) monically excited spherical bubble model. Period. Polytech. 78. Manchein, C., da Silva, R.M., Beims, M.W.: Proliferation of Mech. Eng. 56(2), 125 (2012) stability in phase and parameter spaces of nonlinear systems. 58. Hos, ˝ C., Champneys, A.R.: Grazing bifurcations and chat- Chaos 27(8), 081101 (2017) ter in a pressure relief valve model. Phys. D 241(22), 2068 79. da Silva, R.M., Manchein, C., Beims, M.W.: Controlling (2012) intermediate dynamics in a family of quadratic maps. Chaos 59. Hegedus, ˝ F.: Stable bubble oscillations beyond Blake’s crit- 27(10), 103101 (2017) ical threshold. Ultrasonics 54(4), 1113 (2014) 80. Holzfuss, J., Rüggeberg, M., Mettin, R.: Boosting sonolu- 60. Hegedus, ˝ F., Klapcsik, K.: The effect of high viscosity on minescence. Phys. Rev. Lett. 81, 1961 (1998) the collapse-like chaotic and regular periodic oscillations of 81. Krefting, D., Mettin, R., Lauterborn, W.: Two-frequency a harmonically excited gas bubble. Ultrason. Sonochem. 27, driven single-bubble sonoluminescence. J. Acoust. Soc. Am. 153 (2015) 112(5), 1918 (2002) 61. Hegedus, ˝ F.: Topological analysis of the periodic structures 82. Yasuda, K., Torii, T., Yasui, K., Iida, Y., Tuziuti, T., Naka- in a harmonically driven bubble oscillator near Blake’s criti- mura, M., Asakura, Y.: Enhancement of sonochemical reac- cal threshold: infinite sequence of two-sided Farey ordering tion of terephthalate ion by superposition of ultrasonic fields trees. Phys. Lett. A 380(9–10), 1012 (2016) of various frequencies. Ultrason. Sonochem. 14(6), 699 62. Parlitz, U., Lauterborn, W.: Resonances and torsion numbers (2007) of driven dissipative nonlinear oscillators. Z. Naturforsch. A 83. Merouani, S., Hamdaoui, O., Rezgui, Y., Guemini, M.: Sen- 41(4), 605 (1986) sitivity of free radicals production in acoustically driven bub- 63. Parlitz, U., Lauterborn, W.: Period-doubling cascades and ble to the ultrasonic frequency and nature of dissolved gases. devil’s staircases of the driven van der Pol oscillator. Phys. Ultrason. Sonochem. 22, 41 (2015) Rev. A 36(3), 1428 (1987) 84. Batson, W., Zoueshtiagh, F., Narayanan, R.: Two-frequency 64. Lai, Y.C.: Driving trajectories to a desirable attractor by excitation of single-mode Faraday waves. J. Fluid Mech. using small control. Phys. Lett. A 221(6), 375 (1996) 764, 538 (2015) 65. Gadaleta, S., Dangelmayr, G.: Optimal chaos control 85. Sahoo, B., Panda, L.N., Pohit, G.: Two-frequency parametric through reinforcement learning. Chaos 9(3), 775 (1999) excitation and internal resonance of a moving viscoelastic 66. Martínez-Zérega, B.E., Pisarchik, A.N., Tsimring, L.S.: beam. Nonlinear Dyn. 82(4), 1721 (2015) Using periodic modulation to control coexisting attractors 86. Sahoo, B., Panda, L.N., Pohit, G.: Combination, principal induced by delayed feedback. Phys. Lett. A 318(1–2), 102 parametric and internal resonances of an accelerating beam (2003) under two frequency parametric excitation. Int. J. Nonlinear 67. Poon, L., Grebogi, C.: Controlling complexity. Phys. Rev. Mech. 78(Supplement C), 35 (2016) Lett. 75, 4023 (1995) 87. Huang, S., Chandre, C., Uzer, T.: Bifurcations as dis- 68. Lai, Y.C., Tél, T.: Transient Chaos. Springer, New York sociation mechanism in bichromatically driven diatomic (2010) molecules. J. Chem. Phys. 128(17), 174105 (2008) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Nonlinear Dynamics Springer Journals

Non-feedback technique to directly control multistability in nonlinear oscillators by dual-frequency driving

Free
21 pages
Loading next page...
 
/lp/springer_journal/non-feedback-technique-to-directly-control-multistability-in-nonlinear-ctvR36YFYO
Publisher
Springer Netherlands
Copyright
Copyright © 2018 by The Author(s)
Subject
Engineering; Vibration, Dynamical Systems, Control; Classical Mechanics; Mechanical Engineering; Automotive Engineering
ISSN
0924-090X
eISSN
1573-269X
D.O.I.
10.1007/s11071-018-4358-z
Publisher site
See Article on Publisher Site

Abstract

Nonlinear Dyn https://doi.org/10.1007/s11071-018-4358-z ORIGINAL PAPER Non-feedback technique to directly control multistability in nonlinear oscillators by dual-frequency driving GPU accelerated topological analysis of a bubble in water Ferenc Hegedus ˝ · Werner Lauterborn · Ulrich Parlitz · Robert Mettin Received: 28 December 2017 / Accepted: 12 May 2018 © The Author(s) 2018 Abstract A novel method to control multistability ing the simulations, the control parameters are the two of nonlinear oscillators by applying dual-frequency amplitudes of the acoustic driving at fixed, commensu- driving is presented. The test model is the Keller– rate frequency pairs. The high-resolution bi-parametric Miksis equation describing the oscillation of a bub- scans in the control parameter plane show that a period- ble in a liquid. It is solved by an in-house initial-value 2 attractor can be continuously transformed into a problem solver capable to exploit the high computa- period-3 one (and vice versa) by proper selection of tional resources of professional graphics cards. Dur- the frequency combination and by proper tuning of the driving amplitudes. This phenomenon has opened Electronic supplementary material The online version of a new way to drive the system to a desired, pre-selected this article (https://doi.org/10.1007/s11071-018-4358-z) attractor directly via a non-feedback control technique contains supplementary material, which is available to without the need of the annihilation of other attractors. authorized users. Moreover, the residence in transient chaotic regimes can also be avoided. The results are supplemented with F. Hegedus ˝ Department of Hydrodynamic Systems, Faculty of simulations obtained by the boundary-value problem Mechanical Engineering, Budapest University of solver AUTO, which is capable to compute periodic Technology and Economics, Budapest, Hungary orbits directly regardless of their stability, and trace e-mail: fhegedus@hds.bme.hu them as a function of a control parameter with the W. Lauterborn · R. Mettin pseudo-arclength continuation technique. Cavitation Bubble Dynamics Group, Drittes Physikalisches Institut, Georg-August-Universität Göttingen, Göttingen, Keywords Control of multistability · Dual-frequency Germany e-mail: werner.lauterborn@phys.uni-goettingen.de driving · Bifurcation structure · GPU programming · Keller–Miksis equation · Nonlinear dynamics R. Mettin e-mail: robert.mettin@phys.uni-goettingen.de U. Parlitz ( ) 1 Introduction Research Group Biomedical Physics, Max Planck Institute for Dynamics and Self-Organization, Göttingen, Germany e-mail: ulrich.parlitz@ds.mpg.de One of the common features of nonlinear systems is multistability [1]; that is, two or more attractors can U. Parlitz coexist at a fixed parameter combination. Multistabil- Institute for Nonlinear Dynamics, ity appears in almost any field of science; for instance in Georg-August-Universität Göttingen, Göttingen, Germany biological systems [2,3], hydrodynamics [4], mechan- 123 F. Hegedus ˝ et al. ical engineering [5], chemical reactions [6,7], neuron There are three main control strategies known in the dynamics [8], climate dynamics [9] or social systems literature, namely non-feedback, feedback and stochas- [10] to name a few. Multistability is also found in bub- tic control [1]. Each has its own advantages and disad- ble dynamics and acoustic cavitation [11], thetopicof vantages. The non-feedback techniques are simple and the present study. In applications, a major challenge in easy to use, where the main aim is to kick the system nonlinear dynamics is how to control multistability and to another stable state [31] or annihilate attractors by switch between coexisting states representing different periodic perturbation (modulation) of a parameter or system behaviour. a state variable [32]. Unfortunately, it cannot be guar- The mechanism of the emergence of multiple attrac- anteed that the system settles down onto the desired tors can be rather different. In conservative systems, for attractor. The kick introduces randomness to the con- example, one can obtain an arbitrarily large number of trol, and extremely long transients are also possible coexisting attractors by introducing small dissipation in the presence of transient chaos. In case of periodic [12] (see also the Kolmogorov–Arnold–Moser (KAM) modulation, there is no direct control over which attrac- theory [13]). The number of attractors is roughly tors are annihilated, since usually the attractor with the inversely proportional to the damping factor [14]. The smallest basin is destructed first. With increasing per- appearance of homoclinic tangencies or the coupling turbation magnitude, the annihilation process contin- of systems with a large number of unstable invariant ues with attractors having successively the next smaller manifolds can also lead to a high level of multistabil- basin. Applying a feedback control strategy, the desired ity. In these cases, even an infinite number of attracting attractor can be targeted directly; moreover, it is a reli- periodic orbits called Gavrilov–Shilnikov–Newhouse able technique to control attractors against external sinks can coexist [15,16]. Another important mecha- noise. However, details are necessary about the state nism leading to multistability is the coupling of a large space [33] or even the solution of the attractor itself number of identical systems, where the number of the [34] or its Jacobian. This makes the feedback control attractors scales with the number of oscillators [17]. technique unusable in certain problems. For instance, in This results in a so-called attractor crowding that makes an acoustically trapped bubble, it is hard to obtain state the system extremely sensitive to external noise lead- space information about the oscillating bubble and feed ing to random hopping between many coexisting states. this information back to the resonator. The third kind Several other mechanisms such as delayed feedback of control (stochastic) also performs attractor annihila- [18], parametric forcing [19] or external noise [20] can tion, this time with the addition of external noise [35]. also induce multistability. The drawback is similar as in case of the non-feedback The importance of the control of multistability can technique: the attractors are destructed in the order of be explained with the fact that different stable states the size of their basin of attraction [1]. That is, the direct represent different system performances. A good exam- control is again lost. Attractors with very small basins ple is the large subject of chaos control [21–25], where are usually buried by the noise immediately, without the main aim is to eliminate unpredictable behaviour. the need of increasing the magnitude of the noise to Furthermore, there are various fields where the con- some higher level. trol of multistability is of paramount importance: laser The main aim of the present study is to propose a pos- technology [26], optical communication [27], cardiol- sible non-feedback control technique capable of driv- ogy [28], genetics [29] or ecology [30]. If multistability ing a sinusoidally excited system to a pre-selected peri- is undesired, then basically two approaches are possi- odic attractor directly. It is based on the observation that ble: make the system monostable or preselect and sta- a period-2 and a period-3 attractor can be continuously bilize an attractor against external noise. In case of a transformed into each other by adding a second, sinu- desirable multistability, where the control between dif- soidal component to the excitation. The proper tuning ferent states (operations) on demand is important, the of the excitation amplitudes is necessary; moreover, task is to properly select an attractor the system should the fixed commensurate frequency combination must approach and/or exclude certain undesired attractors be adjusted according to the periods of the attractors from the dynamics. In these cases, the stabilization being transformed. The advantages and disadvantages against external noise can also be crucial. over other, already existing control techniques, the lim- itations and possible extensions to transformation of 123 Nonlinear oscillators by dual-frequency driving a period-n to a period-m attractor in general, and the and the vapour, p = 3166.8Pa. p denotes the pres- V L application to other nonlinear oscillators are also dis- sure in the liquid at the bubble-liquid interface. The cussed. surface tension is σ = 0.072 N/m, and the liquid kine- −4 The test model is the Keller–Miksis equation well matic viscosity is μ = 8.902 Pa s. The gas inside known in the bubble dynamics community. It is a the bubble obeys a simple polytropic relationship second-order nonlinear ordinary differential equation 3γ describing the oscillation of a bubble in a liquid [11]. 2σ R p = P − p + , (4) G ∞ V The parameter space of the present study is relatively R R large due to the dual-frequency driving (amplitudes and frequencies). The approach to obtain a global where the polytropic exponent is γ = 1.4 (adiabatic behaviour), the equilibrium bubble radius is R = picture about the bifurcation structure in the four- dimensional driving parameter space within reason- 10 μm and the static pressure is P = 1 bar. able time is to exploit the high arithmetic processing For numerical purposes, system (1)–(4) is rewritten power of professional graphics cards (GPUs) using an into a dimensionless form by introducing the dimen- in-house initial-value problem (IVP) solver written in sionless variables C++/CUDA C. The parameter scans are supplemented with results obtained by the boundary-value problem τ = t, (5) 2π solver AUTO [36], which is suitable to compute peri- odic orbits directly regardless of their stability, and fol- y = , (6) low their paths with the pseudo-arclength continuation 2π technique. y = R , (7) R ω E 1 2 Mathematical model and the equations are rearranged to minimize the num- ber of coefficients. The final form of the numerical The test model studied (Keller–Miksis equation [11]) model is describing the evolution of the bubble radius reads ˙ ˙ R R 3 ¨ ˙ 1 − R R + 1 − R y ˙ = y , (8) 1 2 c 3c 2 L L KM R R d ( p − p (t )) y ˙ = , (9) L ∞ 2 = 1 + + , (1) D KM c c dt ρ L L L where R(t ) is the time dependent bubble radius; c = where the numerator, N , and the denominator, D , L KM KM 1497.3m/s and ρ = 997.1kg/m are the sound speed are and density of the liquid domain, respectively. Accord- N = (C + C y ) − C (1 + C y ) ing to the general, dual-frequency treatment, the pres- KM 0 1 2 2 9 2 sure far away from the bubble 1 y y 3 2 2 −C − C − 1 − C y 3 4 9 y y 3 2 1 1 p (t ) = P + P sin(ω t ) + P sin(ω t + θ) (2) ∞ ∞ A1 1 A2 2 − (C sin(2πτ ) + C sin(2π C τ + C )) 5 6 11 12 1 + C y − y C cos(2πτ ) consists of a static component, P , and periodic com- ( ) ( 9 2 1 7 ponents with pressure amplitudes P and P , angular A1 A2 + C cos(2π C τ + C )) , (10) 8 11 12 frequencies ω and ω , and with a phase shift θ.The 1 2 and connection between the pressures inside and outside the bubble at its interface can be written as D = y − C y y + C C . (11) KM 1 9 1 2 4 9 2σ R p + p = p + + 4μ , (3) G V L L The coefficients are summarized as follows: R R where the total pressure inside the bubble is the sum of 1 2σ 2π C = P − p + , (12) 0 ∞ V the partial pressures of the non-condensable gas, p , ρ R R ω L E E 1 123 F. Hegedus ˝ et al. 1 − 3γ 2σ 2π tions forward in time and, after the convergence of the C = P − p + , (13) 1 ∞ V ρ c R R ω transient trajectory, record and save the important prop- L L E E 1 erties of the solutions found. Thereby, our main inter- P − p 2π ∞ V C = , (14) est is in periodic solutions and their organization in the ρ R ω L E 1 reduced driving parameter space discussed below. The 2σ 2π first part (2048 T ) of the integration is regarded as ini- C = , (15) ρ R R ω L E E 1 tial transient and dropped followed by integration of 4μ 2π the system further by 8192 T for the proper conver- C = , (16) ρ R ω 1 gence of average quantities like Lyapunov exponent or winding number, where T is the period of the dual- P 2π A1 C = , (17) frequency excitation function. The scheme employed ρ R ω L E 1 is the explicit and adaptive Runge–Kutta–Cash–Karp P 2π A2 method with embedded error estimation of orders 4 C = , (18) ρ R ω L E 1 and 5 (the algorithm is adapted from [38]). The quan- tities stored are the points of the (global) Poincaré sec- ω P 2π 1 A1 C = R , (19) 7 E tion (see the choice in Sect. 4 below), which are stan- ρ c R ω L L E 1 dard ingredients for bifurcation diagrams; furthermore, ω P 2π 1 A2 the maximum bubble radius and the maximum absolute C = R , (20) 8 E ρ c R ω L L E 1 value of the bubble wall velocity, which are important R ω E 1 for applications; finally, the period,the Lyapunov expo- C = , (21) 2π c L nent and the winding number of the attractors found, C = 3γ, (22) quantities that are essential for a detailed analysis of bifurcation structures. C = , (23) A strategy to represent the results of parametric stud- ies involving high-dimensional parameter spaces con- C = θ. (24) sists in creating high-resolution bi-parametric plots, a rapidly spreading technique in the investigation of The angular frequencies ω and ω are normalized 1 2 nonlinear systems with a high-dimensional parame- by the linear, undamped angular eigenfrequency [37] ter space [39–49]. The system studied here, a bub- ble in water with dual-frequency acoustic excita- 3γ( P − p ) 2(3γ − 1)σ ∞ V ω = − (25) tion, has a four-dimensional driving parameter space 2 3 ρ R ρ R L L E E ( P , P ,ω ,ω ). For simplicity, the phase shift is A1 A2 R1 R2 set to θ = 0. It is restricted in the present study by the of the unexcited system that defines the relative fre- requirement that the relative-frequency pairs obtained quencies as from ω and ω are rational. The pressure amplitudes R1 R2 1 P and P are taken as the main control parameters at A1 A2 ω = , (26) R1 ω fixed relative-frequency pairs. The range of each pres- sure amplitude is 0–5 bar, investigated at first with a ω = . (27) R2 resolution of 501 steps for an overview. The relative- frequency values ω and ω are chosen from the R1 R2 With R = 10 µm and with the other constants in (25) following set: given before in this section, the eigenfrequency is f = ω /2π = 340 kHz. 1 1 1 1 1 2 3 5 10 , , , , , , , , . (28) 10 5 3 2 1 1 1 1 1 3 Numerical implementation and parameter choice Exploiting symmetries in the driving parameter space, this gives i = 36 combinations of two fre- i =1 The simplest technique to solve the system (8)–(9)isto quencies. Observe that two orders of magnitude dif- take an initial-value problem solver, integrate the equa- ference in the frequency range are covered and that, 123 Nonlinear oscillators by dual-frequency driving for instance, ω = 1/5 and ω = 1/1 means that R1 R2 T T 1 2 =2 R2 in an experiment, the bubble is driven by 0.2 ω and by the resonance frequency ω . Computations are per- formed at every possible frequency combination of the above set, which means altogether 36 two-dimensional =3 -1 plots for every quantity to be studied (period, Lya- R1 dual-frequency excitation punov exponent, etc). Moreover, owing to the fine -2 0 0.5 1 1.5 2 2.5 3 details in the diagrams, high-resolution surveys are [-] necessary. In order to reveal possible coexisting attrac- tors, 10 randomized initial conditions are used at each Fig. 1 Example for graphs of a dual-frequency excitation (red) and its two harmonic components as mainly used in the study parameter set. Indeed, the 10 different trajectories can (ω = 3, black; ω = 2, blue) with pressure amplitudes R1 R2 converge to different attractors with distinct proper- P = P = 1 bar. The vertical lines marked by T (black A1 A2 1 ties (e.g. period). Altogether, a single plot contains solid), T (blue solid) and T (red dashed) designate the periods 501 × 501 × 10 ≈ 2.5 million initial-value prob- of the high-frequency component, the low-frequency component lems. and the dual-frequency driving, respectively. T = 3 × T = 2 × T = 3. (Color figure online) Even a single bi-parametric plot requires large com- putational capacities; thus, the high processing power of professional videocards (GPGPUs) is exploited. During the parameter studies, millions of simple 4 Global Poincaré section and the period-reducing (only second order) and independent ordinary dif- phenomenon ferential equations are to be solved. Moreover, the algorithm employed needs only function evaluations. The parameters investigated in the present study are Therefore, our problem is well suited for paralleliza- the pressure amplitudes ( P , P ) and the relative fre- A1 A2 tion on GPUs. The program code is written in a quencies (ω ,ω ) of the dual-frequency excitation R1 R2 C++ and CUDA C software environment. The avail- with phase shift θ = 0, i.e. C = 0. The main con- able GPUs are a Titan Black card (Kepler archi- trol parameters are the pressure amplitudes at several tecture, 1707 GFLOPS double-precision processing fixed frequency combinations. The ratio of each fre- power), two Tesla K20 cards (Kepler architecture, 1175 quency pair is rational; therefore, periodicity of the GFLOPS) and eight Tesla M2050 cards (Fermi archi- dual-frequency driving is guaranteed (quasiperiodic tecture, 515 GFLOPS). The application of double- forcing is excluded). precision floating-point arithmetic is mandatory in bub- As an example, a dual-frequency excitation is pre- ble dynamics due to the possibly large-amplitude, sented in Fig. 1 as a function of time together with its collapse-like response of the system [11]. The final, two periodic components at ω = 3 and ω = 2. For R1 R2 optimized code is approximately about 50 times faster definiteness, the pressure amplitudes are set to P = A1 on the Titan Black, about 30 times faster on the Tesla P = 1 bar. The black and blue curves are the high- A2 M2050 and about 120 times faster on the Tesla K20 and low-frequency components, respectively. Since the than on a four-core Intel Core i7-4790 CPU. Sur- dimensionless time τ is defined by means of the first fre- prisingly, the professional Tesla K20 card was more quency component, see Eq. (5), the normalized period than two times faster than the gamer Titan Black card corresponding to the relative frequency ω is T = 1, R1 1 even though the theoretical peak processing power is indicated by the vertical black line in Fig. 1. The nor- higher for the Titan Black GPU. A thorough perfor- malized period of the second component (relative fre- mance analysis between GPUs and CPUs solving large quency ω ) can be calculated with the frequency ratio; R2 numbers of ordinary differential equations is beyond that is, T = 1/C = ω /ω = 3/2 = 1.5 (blue 2 11 R1 R2 the scope of the present study. The interested reader vertical line). Compare these periods with the argu- is referred to the publications [38,50–52]for more ments of the harmonic functions in equation (10) and details. keep in mind that the phase shift is zero (θ = C = 0). As only normalized periods will be used throughout the article, this very often occurring quantity will be simply called “period”. p , p [bar] A1 A2 F. Hegedus ˝ et al. The sum of the two harmonic components results is no restriction to small-amplitude perturbations cor- in a still periodic, but non-harmonic forcing function responding to the second sinusoidal component. presented by the red curve in Fig. 1. The period of this function is T = 3T = 2T = 3, whichisusedasa 5.1 The branching phenomenon 1 2 global Poincaré section for the system (8)–(9). Con- sequently, during the numerical integration of the sys- Out of the data sets for the 36 relative-frequency com- tem with an initial-value problem solver, the continu- binations, four diagram pairs are shown in Fig. 2 at the ous trajectories are sampled at time instants τ = n · T four different relative-frequency pairs (ω ,ω ) = n R1 R2 (n = 0, 1, 2,...). Accordingly, the period, T ,ofa (0.2, 0.1), (1, 0.1), (3, 1) and (3, 2) organized in rows periodic orbit is defined as T = N · T , N ∈ IN, and and for two different quantities (Lyapunov exponent the solution is called period- N orbit. Observe, however, and period) organized in columns. The left column that in the special cases of P = 0 bar or P = 0 bar, contains the Lyapunov-exponent diagrams, where the A1 A2 the period T , which defines the global Poincaré section, colour-coded area means chaotic solutions (positive does not coincide with the period of the actual driving exponent) and where in the greyscale domains there T or T .For P = 0 bar, the dual-frequency forc- are periodic attractors (negative exponent). In the right 2 1 A1 ing becomes a purely harmonic function shown by the column, the periods of the (converged) solutions are blue curve in Fig. 1, and the continuous trajectories presented up to period-6. Periodic solutions higher than are sampled only after every second period of excita- six including chaos can be found in the black domains. tion, 2T . Therefore, an originally period-2k solution In the upper two rows, the plots correspond to relative is observed only as a period-k orbit. Similarly, when frequencies lower than or equal to the linear resonance P = 0 bar, the trajectories are sampled only after frequency (ω ≤ 1), while the lower two rows have A2 R1,2 every third period of the excitation, 3T (compare the relative frequencies higher than or equal to the linear black and red vertical lines in Fig. 1), and the period resonance frequency (ω ≥ 1). R1,2 of an originally period-3k orbit is reduced from 3k to Since detailed investigations of dual-frequency driven k.This period-reducing phenomenon, occurring when systems for arbitrary driving amplitudes with ratio- one of the two components of the forcing function is nal frequency ratio are absent in the literature, only vanishing and only a purely harmonic excitation is left, very few information exists about the bifurcation struc- plays an important role in the bifurcation structure of ture shown in Fig. 2. It can be clearly seen that the the system, discussed in more detail in Sect. 5; further- dual-frequency driving causes a very complex interplay more, period reduction necessarily takes place for other between the resonances originating from the single- frequency pairs as long as the ratio of a pair is rational. frequency driven system (horizontal and vertical axes in Fig. 2), see especially the first two rows of the figure. 5 Results from a global scan Their comprehensive topological analysis is beyond the scope of the present study; instead, it intends to add Throughout the following subsections, the fundamental some pieces to our knowledge by looking for special mechanism capable of controlling multistability in har- features in the diagrams. In particular, we focus on a monically driven nonlinear oscillators is presented. The very specific phenomenon, the branching mechanism, basis is the period-reducing phenomenon discussed turning out to be a special feature of dual-frequency above, by which the dual-frequency treatment can gen- driven nonlinear oscillators with rational frequency erate “bridges” between different periodic attractors. combinations, which can be used for the control of mul- First, a global picture is shown via a series of high- tistability. resolution bi-parametric plots of Lyapunov-exponent Observe that on the horizontal axis in the third row and periodicity diagrams at different frequency com- and second column of Fig. 2, three period-3 branches binations (the control parameters are the excitation emerge from the segment marked by P3B3 (period: 3, amplitudes). Next, the details are demonstrated through number of the branches: 3). Similarly, three branches a specific example with frequencies ω = 3 and appear from the segments P4 B3 and P3B3 on the hor- R1 ω = 2 (see again Fig. 1 for the driving signals). izontal axis, and two branches are merged together at R2 It must be emphasized that the amplitudes for both fre- segments P2B2 and P3B2 on the vertical axis in the quencies of the driving can be arbitrary; that is, there last row and second column of Fig. 2, better to be seen 123 Nonlinear oscillators by dual-frequency driving Fig. 2 List of high-resolution bi-parametric plots as a function of the pressure amplitudes at different relative-frequency pairs organized in rows. Left column: Lyapunov- exponent diagrams separating the chaotic and periodic solutions. Right column: Periodicity diagrams, where periodic solutions are plotted with different colours up to period-6. Period-7 or higher solutions (including chaos) can be found in the black regions. The pressure amplitudes P and P are A1 A2 given in bar in the enlargement and extension in Fig. 3. Observe In order to get a deeper insight into the branch- also that the number of the branches is exactly the ing mechanism at the relative-frequency pair ω = R1 same as the value of the corresponding relative fre- 3,ω = 2, computations were repeated with much R2 quency along the respective axis. This phenomenon higher resolution of the pressure amplitudes (Fig. 3). will be referred to as branching mechanism through- The ranges of the control parameters are P ∈ A1 out the paper and will be investigated in more detail for [0, 8] (bar) and P ∈[0, 5] (bar) with resolutions A2 the case ω = 3,ω = 2. 4001 and 2501, respectively. The number of initial R1 R2 123 F. Hegedus ˝ et al. Fig. 3 Bi-parametric periodicity diagram at the relative- grid in the y – y phase plane. The colour code is the same as in 1 2 frequency pair ω = 3,ω = 2 with increased resolution the right column of Fig. 2. The pressure amplitudes P and P R1 R2 A1 A2 (4001 × 2501) of the control parameters P and P . The num- are given in bar. (Color figure online) A1 A2 ber of the initial conditions is 5×5 = 25 defined on an equidistant conditions at a control parameter pair ( P , P ) is bordered each by a yellow (period-4) band at the right A1 A2 increased from 10 to 5 × 5 = 25 defined on an equidis- border (a period-doubled zone) followed by a thin black tant grid on the y – y phase plane. Thus, 2501×4001× line (indicating period-8 and higher periods including 1 2 25 ≈ 250 million initial-value problems have to be chaos that are encoded that way). Similarly, the three solved. This huge amount of computations is divided dark blue P3 bands are bordered each by a light blue into 5 × 8 = 40 (Δ P = 1 bar) segments and dis- (period-6) band, again indicating the start of a period- A1,2 tributed to different GPUs. doubling sequence. In addition, the three yellow P4 At first sight, Fig. 3 contains only slightly more bands are bordered by a thin black line, presumably information about the branching mechanism than the also the beginning of a period-doubling sequence to be bottom right diagram in Fig. 2, since it is hard to rep- seen only at higher resolution. Note that also other dark resent the many coexisting attractors in a single plot blue areas are bordered on one side by a period-doubled obtained by the increased number of initial conditions zone of light blue followed by black. Thus, the branches (in case of coexistence, the solution with the high- themselves are period-doubling objects. It should also est period is plotted). Investigating the system period- be noted that in any case, two of the three branches in by-period, however, this high-resolution bi-parametric each Pn B3, n = 2, 3, 4, turn to the left (lower driving scan becomes very helpful to understand the branching amplitudes P ) and one to the right (higher driving A1 mechanism (explored rigorously in the next sections); amplitudes P ). A1 in particular, to identify the connected branches with different periods confidently, such as the ones originat- ing from P4B3, P3B3, P2B3, P3B2 or P2B2, where, 5.2 Coexisting period-1 solutions for instance, P2B3 and P2B2 are connected. Typical colour combinations can be seen at the bor- Although the red period-1 domain in Fig. 3 shows no der of the branch families. The three green P2 bands are sign of multiple branches, the examination of its inner 123 Nonlinear oscillators by dual-frequency driving structure—the organization of the coexisting period- 1 orbits—reveals important features of the branching mechanism. They turn out to be general for other peri- ods as well. By filtering the results for period-1 orbits and by counting the number of the different attractors at a given control parameter pair, the structure of the coexisting period-1 solutions can be made visible in the P – P plane, see Fig. 4 with a magnification of A1 A2 the zone near the origin in the bottom panel. Domains with different colours represent areas with a differ- ent number of coexisting attractors up to four. These domains are bounded by saddle-node (SN, solid) and period-doubling (PD, dashed) curves computed by the boundary-value problem solver AUTO introduced in more detail later. The boundaries computed by AUTO are continuous curves. Some fractal-like shapes in the colour-coded domains are an artefact of the lack of a sufficient number of initial conditions or the presence of extremely long transients during the initial-value- problem computation Crossing one of the boundaries means the increment or decrement of the number of the attractors by one. Fig. 4 Number of coexisting period-1 orbits in the P – P A1 A2 Multiple stable period-1 orbits in a bi-parametric parameter plane at the relative-frequency pair ω = 3,ω = R1 R2 space usually originate via cusp bifurcations forming 2. The saddle-node (SN, solid) and period-doubling (PD, dashed) pairs of SN points (hysteresis) and giving rise to hys- curves, computed by the boundary-value problem solver AUTO, are the boundaries of the domains with a different number of teresis zones. The overlapping of such hysteresis zones coexisting attractors. The pressure amplitudes P and P are A1 A2 may produce domains with a high number of coexist- giveninbar ing period-1 attractors [53]. The sole cusp bifurcation in the bottom panel of Fig. 4, however, cannot explain the mosaic-type pattern of the period-1 structure. The cation points are indicated by PF, SN and PD, respec- investigation of the two special cases, where one of the tively. amplitudes of the dual-frequency driving P or P The period-reducing phenomenon discussed in detail A1 A2 is zero (vertical or horizontal axis in Fig. 4), can help in Sect. 4 plays an important role in the period-1 bifur- to understand this pattern. cation structure. Observe that the applied frequency combination is the same as in Fig. 1; furthermore, the first component of the dual-frequency driving is zero. Consequently, the driving is a purely harmonic 5.3 The limiting case P = 0 and the function proportional to the blue curve in Fig. 1, and A1 period-reducing phenomenon an originally period-2k solution is observed only as period-k orbit. Thus, the first bifurcation of the sta- The bifurcation structure with P = 0 bar is shown in ble period-1 curve originating from the equilibrium A1 Fig. 5, where the first component of the Poincaré sec- solution y (both amplitudes are zero) is a pitchfork 1 E tion Π( y ) is plotted as a function of P . The colour point PF, instead of a period-doubling point PD, which 1 A2 code is the same as in the cases of Fig. 3 and the right gives birth of two period-1 branches 2 × P1. Only the column of Fig. 2. For instance, red, green and blue subsequent bifurcations exhibit a Feigenbaum cascade curves are period-1, -2 and -3 orbits, respectively. Some (2 × P2 → 2 × P4 → ··· ). Similarly, the first bifur- of the periodic orbits are also marked by N ×PX, where cation point in a cascade must be a symmetry-breaking N is the number of the coexisting attractors of period X. bifurcation in case of a symmetric orbit in a symmetric The pitchfork, saddle-node and period-doubling bifur- nonlinear oscillator [54]. 123 F. Hegedus ˝ et al. Fig. 5 One-dimensional bifurcation structure with P = 0at periodic orbits are marked by N × PX,where N is the number of A1 the relative-frequency pair ω = 3, ω = 2, where the first the coexisting attractors of period X. The pitchfork, saddle-node R1 R2 component of the Poincaré section Π( y ) is plotted as a function and period-doubling bifurcation points are marked by PF, SN of the control parameter P . The colour code is the same as in and PD, respectively A2 the cases of Fig. 3 and the right column of Fig. 2. The relevant 1.4 Two examples of bubble radius versus time curves 1.3 y (τ ) are plotted in the top panel of Fig. 6 so as to 1.2 help understanding the period-reducing mechanism. The blue and black period-1 solutions are initiated 1.1 at P = 1 bar from the black dots marked by the A2 bold numbers 1 and 2 in Fig. 5, respectively. The blue 0.9 harmonic function in the bottom panel in Fig. 6 is 0.8 the purely harmonic driving of the system, see also 0.7 Fig. 1. It is clear that if the Poincaré sections are taken at time instances τ = n · T (n = 0, 1, 2,...), 0.6 0 0.5 1 1.5 2 2.5 3 which is used during the present simulations, the solu- [-] tions are regarded as period-1 orbits. In contrast, if the =2 R2 Poincaré sections are taken at time instances τ = n · T 1 n 2 (n = 0, 1, 2,...), which is a common choice in case of a purely harmonic driving, the solutions are period-2 -1 orbits. To prevent ambiguity in the definition of the peri- odicities, solutions will be referred to as “originally” -2 0 0.5 1 1.5 2 2.5 3 period- N orbits if the latter definition of the Poincaré [-] section is applied. Although the blue and black attrac- tors in Fig. 6 are different (each has its own basin of Fig. 6 The period-reducing phenomenon. Top panel: Bubble radius versus time curves y (τ ) (pressure amplitude P = 0bar 1 A1 attraction), the only difference between them is a phase P = 1 bar) starting at the black dots at τ = 0marked by A2 shift in time by one period of the driving (Δτ = 1.5). the bold numbers 1 and 2 in Fig. 5. Bottom panel: The sig- Altogether two types of bifurcation scenarios of nal of the single-frequency driving with period τ = T .The periodic orbits are observable in Fig. 5. For odd basic global Poincaré section applied during the present computations is defined by T (dashed, red vertical line). (Color figure online) periods (period before a bifurcation cascade takes place), the first point is a pitchfork bifurcation followed by a Feigenbaum period-doubling cascade (PF → y [-] p [bar] A2 Nonlinear oscillators by dual-frequency driving (A) (B) 2 0.7 =3 R1 0.6 =2 R2 1 1.5 P1 2 0.5 P =0 bar E 1 A1 PD P1 +P1 1 2 0.4 P =0.3 bar 0.3 E,u A1 E SN P =0.6 bar P1 P1 A1 0.2 0.5 PD P1 PF 0.1 E,u 2 P1 +P1 P =0 bar 1 2 A1 -0.1 2E -0.5 -0.2 012345 0.15 0.2 0.25 P [bar] P [bar] A2 A2 Fig. 7 Period-1 bifurcation curves computed by the boundary- node (SN) bifurcation points, respectively. The left panel (A) value problem solver AUTO. The black and red curves are stable shows the case with P = 0, while the right panel (B) A1 and unstable solutions, respectively. The asterisk, crosses and presents results at pressure amplitude P = 0.3 bar and at A1 dots are the pitchfork (PF), period-doubling (PD) and saddle- P = 0.6 bar. (Color figure online) A1 PD → PD → ··· ). Some examples are the already In order to illustrate this effect, it is useful to solve investigated period-1 curve P1 starting from y ,the and compute complete bifurcation curves of periodic 1 E period-3 orbits P3 originating from a SN point at orbits with a boundary-value-problem solver. As an P = 1.5 bar or the period-5 branches P5 appearing example, in Fig. 7A, the second component of the A2 again via a SN point at P = 2 bar. Observe that the Poincaré section Π( y ) of the period-1 curves is plot- A2 2 PF points have no influence to the colour, i.e. period. ted as a function of the control parameter P com- A2 In case of even basic periods, the period is immedi- puted by the bifurcation analysis and continuation soft- ately halved (period-reducing phenomenon), and the ware AUTO [36]. Here, the first component of the bifurcation scenario is a classic Feigenbaum cascade pressure amplitude P is still zero. AUTO is capable A1 without an initial PF point (PD → PD → ··· ). An of tracking down whole bifurcation curves including example here is the originally period-6 solution 2 × P3 the unstable part even if they contain multiple turn- generated at P = 2.5 bar. ing points, and it can detect the bifurcations and their A2 Through the discussion of Figs. 5 and 6, it has been types. This is the reason why AUTO is commonly used shown that the definition of the Poincaré section has an to study the bifurcation structure of nonlinear systems important influence on the apparent periodicity of the [5,53,55–61]. In Fig. 7, the stable and unstable parts of solutions. However, problems only appear in the two the period-1 branches are indicated by the black and limiting cases of either P = 0 bar or P = 0 bar. In red curves, respectively. The curve originating from A1 A2 the first special case of P = 0 bar, a suitable choice y becomes unstable via a pitchfork bifurcation PF A1 2 E can be either T or T , see the bottom panel of Fig. 6. at P = 0.19 bar producing two stable coexisting 2 A2 Since the addition of any small value to the first pres- period-1 branches. These stable curves lose their stabil- sure amplitude P will destroy the special structure ity approximately at P = 3 bar by a period-doubling A1 A2 introduced in Fig. 5, T will be no longer a period of point PD. the dual-frequency driving but only T . Therefore, in order to keep the definition unified, the usage of T as 5.4 Perturbation of the limiting case P = 0 A1 a global Poincaré section is advantageous. This, at first sight only technical, extension of the definition to the The breakup of the structurally unstable pitchfork bifur- limiting, purely harmonic-driving case is of great help cation due to the effect of a nonzero value of the ampli- in ordering and explaining the bifurcation scenarios of tude of the first harmonic component of the excitation dual-frequency driven systems. ( P = 0.3 bar) is clearly seen in Fig. 7B. The curves A1 E 1 marked by P1 and P1 in Fig. 7A form a single curve 1 2 (y ) [-] (y ) [-] 2 F. Hegedus ˝ et al. Fig. 8 One-dimensional bifurcation structure with P = 0at periodic orbits are marked by N × PX,where N is the number A2 the relative-frequency pair ω = 3, ω = 2, where the first of the coexisting attractors of period X. The saddle-node and R1 R2 component of the Poincaré section Π( y ) is plotted as a function period-doubling bifurcation points are marked by SN and PD, of the control parameter P . The colour code is the same as in respectively. (Color figure online) A1 the cases of Fig. 3 and the right column of Fig. 2. The relevant E ,u E 1 P1 + P1 in Fig. 7B. Likewise, the curves P1 and lower periodicities, and how these multiple orbits are 1 2 1 P1 presented in Fig. 7A form again a single curve separated under the perturbation of the amplitude of the E ,u other harmonic component. In order to understand the P1 + P1 in Fig. 7B composed by a stable and an 1 2 global picture, however, the investigation of the other unstable part separated by a saddle-node bifurcation special case, where the control parameter is P and (black dot) at the turning point. Elevating P from A1 A1 the amplitude of the second component P is zero, 0.3 bar to 0.6 bar, a pair of SN points appears along the A2 E 1 is necessary. The corresponding bifurcation diagram branch P1 + P1 forming a hysteresis. Observe that 1 2 E ,u E 1 2 is shown in Fig. 8, where the first component of the both curves P1 + P1 and P1 + P1 have their 1 2 1 2 Poincaré section Π( y ) is plotted as a function of the own separate evolution as the amplitude P changes; A1 pressure amplitude P . The colour code and the nota- A1 consequently, the black and blue curves in the top of tion system are the same as in Fig. 5.Atfirstsight—in Fig. 6 must become different (not only by a phase shift view of (2) with θ = 0—both cases ( P = 0 and A1 in time) for P > 0 bar as they have their own sepa- A1 P = 0) seem exchangeable. However, switching to A1 rate evolution as well. All the pitchfork structures in the the other frequency component, different bifurcation bifurcation scenarios initiated from odd basic periods scenarios develop (if the frequencies are distinct). are broken up in the same way as presented in Fig. 7B, There are two types of bifurcation scenarios of peri- except for the formation of the hysteresis which not odic orbits. If the basic (original) period is divisible necessarily happens. by 3, then the period reducing takes place immediately decreasing the period of the orbits from 3k to k followed by a Feigenbaum cascade as usual (PD → PD → ··· ). 5.5 The limiting case P = 0 A2 Notable examples are the branches 3 × P1 and 3 × P3 both appearing via SN bifurcations at P = 1.2 bar A1 Shortly in this section, the complex development of and at P = 4.1 bar, respectively. For other basic A1 the period-1 structure in a wide range of the amplitudes periods, there is no period-reducing phenomenon, but P and P will be discussed in detail. Up to now, it is A1 A2 only a period-doubling cascade. Examples of this sce- important only to understand how the period-reducing nario are the period-1 curve departing from the equi- mechanism works for single-frequency driving, how it librium solution y , the period-4 orbits generated at 1 E decomposes periodic solutions into multiple orbits of 123 Nonlinear oscillators by dual-frequency driving 1.4 stable up to the pressure amplitude P ≈ 7 bar. It A1 1.3 loses its stability there via a period-doubling bifur- cation. The originally period-3 family of solutions 1.2 3 × P1 appearing through a SN point is decoupled into 1.1 three period-1 branches existing between the ampli- tudes P ≈ 1.19 bar and P ≈ 4.67 bar. Therefore, A1 A1 0.9 there is a wide range of amplitudes P where four A1 period-1 attractors exist. Figure 9 shows these coexist- 0.8 ing orbits at P = 2 bar starting from the black dots A1 0.7 in Fig. 8. The dashed curve is the period-1 orbit related 0.6 to the curve originating from y . The black, red and 1 E 0 0.5 1 1.5 2 2.5 3 [-] blue solid curves are the orbits of the period-reduced 3 × P1 branches, where the only difference between Fig. 9 The period-reducing phenomenon. Dimensionless bub- them is again a phase shift in time, in this case of one ble radius versus time curves of four coexisting stable period-1 or two periods T (Δτ = 1or2). orbits at pressure amplitude P = 2 bar starting from the black A1 1 dots shown in Fig. 8 at τ = 0. Between the black, red and blue solid curves there is only a phase shift in time with Δτ = 1 or 2. The dashed curve is a degenerate (i.e. repeating only after T instead of after T ), small-bubble-radius-amplitude case that 1 5.6 Perturbation of the limiting case P = 0 A2 exists stably up to P ≈ 7 bar where the oscillation period A1 doubles. (Color figure online) Again, the boundary-value-problem solver AUTO can help to understand the behaviour of the system under a small perturbation of the amplitude of the second, P = 1.5 bar or the period-5 branches initiated at harmonic component of the driving ( P > 0). The A1 A2 P = 4.6 bar. Observe that there are no pitchfork or bifurcation curves related to the four period-1 orbits A1 other special bifurcations present in these cases. can be seen in Fig. 10Afor P = 0 bar. The black A2 Figure 8 explains how the coexistence is possible of and red curves mean stable and unstable orbits, respec- altogether four period-1 stable solutions highlighted in tively. The bifurcation of the saddle-node (SN) and Fig. 4 by the yellow domain. The period-1 orbit P1 period-doubling (PD) points are also detected. They emerges from the equilibrium point y and remains are marked by dots (SN) and crosses (PD). Observe 1 E 2.5 1.4 (B) =3 (A) P =0.1 bar R1 A2 1.2 PD =2 2 1 R2 P1 3 1 P =0 bar A2 1.5 0.8 0.6 SN 0.4 P1 0.5 0.2 P1 3 0 2 -0.2 P1 2E -0.5 -0.4 01234567 8 0.5 1 1.5 2 P [bar] P [bar] A1 A1 Fig. 10 Period-1 bifurcation curves computed by the boundary- while the right panel (B) presents results of the 3 × P1curves value problem solver AUTO. The black and red curves are stable at the pressure amplitude of P = 0.1 bar. The arrows in panel A2 and unstable solutions, respectively. The crosses and dots are the (B) show the motion of the SN bifurcation points, two of them period-doubling (PD) and saddle-node (SN) bifurcation points, towards lower P , one of them towards larger P . (Color figure A1 A1 respectively. The left panel (A) shows the case with P = 0, online) A2 (y ) [-] y [-] 2 1 (y ) [-] 2 F. Hegedus ˝ et al. that the P -ranges of the stable (black) segments in A1 Fig. 10A coincide with the domain of existence of the corresponding parts of the red curves in Fig. 8.Onthe right hand side of Fig. 10, only the evolution of the curves 3 × P1 is presented by increasing the ampli- tude P fromzeroto0.1 bar. It is clear that while the A2 saddle-node bifurcations (black dots) take place at the same value of P when P = 0 bar (see the vertical A1 A2 thin line in Fig. 10B), the locations of these bifurcation points differ when P > 0. This means that each of A2 the three bifurcation curves have their own evolution as the amplitude P changes. Consequently, the black, A2 red and blue solid curves in Fig. 9 become different (not Fig. 11 3D representation of the stable period-1 orbits (Poincaré only by a phase shift in time) for P > 0. A2 section points of the bubble wall velocity) above the parameter plane P – P by filtering the results presented in Fig. 3 by A1 A2 period, here period-1. The period-1 orbits can be separated into four different layers marked by the blue, black, red and green surfaces. The indication of the curves in the Π( y )– P (Fig. 10) 2 A1 6 Period 1: the complete picture m or Π( y )– P (Fig. 7)planesis P1 ,where n is the original 2 A2 period and m = 1, ..., n is a serial number, P1 means period- 1 solution. The black surface originating from the equilibrium Due to the exhaustive discussion in the preceding sec- E 1 solution of the system is marked by P1 (instead of P1 to 1 1 tions, it is clear now how the system behaves under show the connection with the equilibrium point). (Color figure single-frequency driving when a second, harmonic online) component is added as a perturbation. The huge amount of results presented in Fig. 3 can help to establish a global overview of the evolution of all the period-1 Figure 11 indicates a very interesting phenomenon, solutions by means of filtering the results by period, namely the originally period-2 solutions on P1 can here period-1. Figure 11 shows these period-1 solu- be transformed into the originally period-3 solutions tions represented in a three-dimensional plot, where the on P1 through the blue surface. Similarly, orbits on 2 2 second component of the Poincaré section, Π( y ),the P1 can be transformed into orbits on P1 via the 2 3 bubble wall velocity, is plotted as a function of the two green surface. These are smooth transformations by main control parameters P and P . The different means of the continuous change of the amplitudes P A1 A2 A1 period-1 orbits can be decomposed into several layers and P of the dual-frequency driving. Naturally, all A2 indicated by the blue, black, red and green surfaces. these orbits are considered as period-1 attractors due In the limit cases when one of the pressure amplitudes to the specific choice of the global Poincaré section. is zero (planes Π( y )– P or Π( y )– P ), the points Consider, however, the following scenario of parame- 2 A1 2 A2 are plotted with bigger markers to be able to distin- ter variations visualized by the supplementary movie guish such special solutions more easily. These solu- file (Online Resource 1) where the top panel shows the tions, which form complete bifurcations curves pre- dimensionless bubble radius curves y as a function of sented in Figs. 7A and 10A, are marked by P1 , where the dimensionless time τ , and the bottom panel repre- n is the original period and m = 1,..., n is a serial sents the dual-frequency driving signal P (τ ) (exclud- number, P1 means period-1 solution. When both P ing the static ambient pressure). A1 and P are zero then the system is in equilibrium and Firstly, let us start the investigation with a solu- A2 the dimensionless bubble wall velocity and its Poincaré tion lying on the bifurcation curve 2 × P1inFig. 5 1 2 section y are zero (it is the origin of the diagram). ( P1 or P1 in Fig. 11) so that the amplitude of 2 E 2 2 The black period-1 surface originating from this point is the single-frequency excitation is somewhere between E 1 marked by P1 (instead of P1 , emphasizing the equi- P = 1.4 bar and P = 2.9 bar ( P = 0 bar) with A2 A2 A1 1 1 librium origin E). There is only one curve of this type, relative frequency ω = 2. Specifically, in Online R2 which is a surface of small-amplitude bubble oscilla- Resource 1, the pressure amplitude is chosen to be tions. P = 2.5 bar. In this video, besides the originally A2 123 Nonlinear oscillators by dual-frequency driving Fig. 12 A 3D representation of the bifurcation structure, com- limit cases of P = 0bar or P = 0 bar. Panel B is an enlarge- A1 A2 puted by AUTO, of the blue and black layers presented in Fig. 11. ment of the near-origin zone to better show the bifurcation-point The black and red curves are the stable and unstable solutions, area. (Color figure online) respectively. The blue and green curves form the skeleton of the period-2 (red) attractor, which is the initial state of the tinuous domain in the pressure amplitude—frequency transformation, an originally period-3 (blue) attractor parameter plane [53,61]. Therefore, a smooth transfor- coexists regarded as a final target of the control prob- mation of a period-2 solution into a period-3 orbit of a lem. The colour code of these attractors is the same single-frequency driven system is possible by a tempo- both in Fig. 5 and in the top panel of Online Resource rary addition of a second harmonic component to the 1; in addition, by comparing the top and bottom panels, driving. This is an efficient way to control multistability the periodicities of the attractors are clearly visible. in harmonically driven nonlinear oscillators. Accord- Secondly, guide the period-2 orbit smoothly to the ing to the best of the authors’ knowledge, this phe- 1 2 curve 3 × P1 shown in Fig. 8 (either P1 or P1 in nomenon has not yet been published in the literature. 3 3 Fig. 11) by varying the pressure amplitudes continu- During this third step in Online Resource 1, the change ously. By the end of the second operation, the sys- of P is not necessary and kept constant. Observe how A1 tem is again driven by a single frequency (ω = 3) the initial and the final parameter set become identi- R1 with pressure amplitude between P = 1.2 bar and cal; that is, single-frequency driving with amplitude A1 P = 4.7 bar ( P = 0 bar). In Online Resource 1, P = 2.5 bar at ω = 2. The only difference is A1 A2 A1,2 R1,2 this transformation takes place along a circle in the the state of the system: operation on the period-3 (blue) P – P plane with radius 2.5 bar; that is, the final instead on the period-2 (red) attractor. A1 A2 value of the first pressure amplitude is P = 2.5 bar. In order for a better visualization of the surfaces pre- A1 The continuously changing trajectory during the trans- sented in Fig. 11, hundreds of bifurcation curves are formation is marked by the black curve in the top panel. computed by AUTO departing from the results shown Thirdly, keep the single-frequency driving but change in Fig. 10A. These surfaces, composed by lines, are both the amplitude P and the relative frequency ω summarized in Figs. 12, 13 and 14. The stable and A1 R1 of the driving (in general) to transform the solution unstable parts are the black and red curves, respectively. from the red curve 3 × P1inFig. 8 to the blue curve As a skeleton, all the bifurcation curves of the limit P3inFig. 5. This third step is possible, since both the cases ( P = 0 bar or P = 0 bar) are presented by A1 A2 3 × P1 and the P3 curves are related to the same sub- the blue (stable) and the green (unstable) curves, which harmonic resonance of order 1/3, which forms a con- can help to identify the surfaces more easily (compare 123 F. Hegedus ˝ et al. (A) (B) 2.5 1.5 1.5 1 0.5 0.5 -0.5 0 3 SN P1 5 PF 4 P1 8 8 1 5 1 4 0 0 P [bar] 0 P [bar] 1 A2 A2 P [bar] P [bar] A1 A1 Fig. 13 A 3D representation of the bifurcation structure, com- The blue and green curves form the skeleton of the limit cases puted by AUTO, of the red layer presented in Fig. 11. The black of P = 0bar or P = 0 bar. (Color figure online) A1 A2 and red curves are the stable and unstable solutions, respectively. (A) (B) 2.5 1.5 1.5 0.5 0.5 -0.5 P1 2 2 P1 4 -0.5 P1 P1 3 3 3 PF 3 SN 2 P1 0 0 0 P [bar] 0 P [bar] A2 P [bar] P [bar] A1 A1 A2 Fig. 14 A 3D representation of the bifurcation structure, com- The blue and green curves form the skeleton of the limit cases puted by AUTO, of the green layer presented in Fig. 11. The black of P = 0bar or P = 0 bar. (Color figure online) A1 A2 and red curves are the stable and unstable solutions, respectively. also with Figs. 7A, 10A). As usual, the pitchfork (PF), the right hand side is a magnification of the interesting saddle-node (SN) and period-doubling (PD) bifurca- parts. tion points are marked by asterisk, dots and crosses, From Fig. 12, it is clear that the blue and black sur- respectively. The left hand side of these figures repre- faces in Fig. 11 are connected for very small values sents the surfaces in the full parameter domain, while of P . The separation takes place approximately at A1 P = 0.33 bar via a cusp bifurcation forming a hys- A1 (y ) [-] (y ) [-] 2 2 (y ) [-] (y ) [-] 2 Nonlinear oscillators by dual-frequency driving panel of Fig. 15, there are two branches along which the following conversion occurs: 2 × P3 = P6 ↔ 3 × P3 = P9; that is, an originally period-6 solution is transformed into an originally period-9 orbit. Simi- larly, the bottom panel of Fig. 15 shows four bands con- verting originally period-8 orbits into period-12 orbits (2 × P4 = P8 ↔ 3 × P4 = P12). As a generaliza- tion, the present study found that applying the relative- frequency combination ω = 3 and ω = 2, a trans- R1 R2 formation between period-2k and period-3k attractors is possible, here k = 1, 2,.... 8 Winding numbers and compatibility of periodic orbits Naturally, not every kind of period-2k and period-3k orbits can be transformed into each other; in other words, there is a “compatibility” condition. Observe that the first PD points of the period-1 branches in Fig. 5 (2 × P1) and in Fig. 8 (3 × P1) are connected via the period-doubling curve PD shown in the top panel of Fig. 4. It is well known in the literature that the quan- tity called winding number w = nk/mk is invariant and Fig. 15 Number of the coexisting period-3 orbits (top panel) and period-4 orbits (bottom panel) in the P – P parameter plane A1 A2 rational near a period-doubling or a saddle-node point at relative frequencies ω = 3and ω = 2 R1 R2 [62,63]. Here, mk is the period of a given attractor and nk is the number of the rotations of a nearby trajectory around the attractor during one period. Therefore, to teresis, see also the bottom panel of Fig. 4. An example each bi-parametric PD or SN bifurcation curve a sin- for this hysteresis, which exists between the regions of gle, rational winding number can be associated. When a P > 0.33 bar and P < 1.19 bar, is already intro- A1 A1 transformation between period-2k and period-3k orbits duced in Fig. 7Bfor P = 0.6 bar. As expected, the A1 is possible, then it must be possible between their cor- red, standalone surface in Fig. 11 has no connections responding bifurcation points as well, see, for instance, with other surfaces, see Fig. 13. The green surface in the PD curve in Fig. 4 which has the winding number Fig. 11 together with its unstable counterpart forms a w = 1/2. Consequently, the transformation is possible simple banded, dual-layer surface presented in Fig. 14. only if the winding numbers of such bifurcation points are compatible; that is, when they are equal. 7 Higher periods 9 Discussion If the smooth transformation works between certain period-2 and period-3 orbits, there must be other com- According to the literature, the only option to directly binations of periodic solutions where the transforma- target an attractor in case of multistability is to apply a tion is possible. In order to investigate these other pos- feedback control technique. The main drawback of the sibilities, the number of the coexisting period-3 and feedback control methods is that explicit knowledge of period-4 attractors is plotted in the top and the bottom the state space (modified targeting method [33], bush- panels of Fig. 15, respectively. The colour codes are the like paths to a pre-selected attractor [64], reinforcement same as in case of Fig. 4 (period-1 orbits). Let us focus learning [65]) or the properties of the targeted attrac- only on those parts of the complex bifurcation struc- tor (trajectory selection by a periodic feedback [34]) is tures where the transformations take place. In the top necessary. Control of multistability is possible by delay 123 F. Hegedus ˝ et al. feedback as well, which can be combined with chaos cation structure of many nonlinear oscillators have the control easily [66]. In that case, the feedback stabi- same basic building blocks [70]. If there is no informa- lizes one of the unstable orbits embedded in a chaotic tion about a state variable, the orientation in the topol- attractor. Unfortunately, it just stabilizes an unstable ogy can be supported; for instance, by measuring and orbit and does not target directly an already existing Fourier transforming an indirect quantity connected to attractor. Parenthetically, the feedback control can sta- the dynamics (in case of bubbles their emitted pres- bilize an attractor against external noise by employing sure signal). The peaks presented in the spectra can the Jacobian of the attractor as a feedback [67]. How- provide valuable information about the nature of the ever, the state of the system must be already near to actual attractor [71], see also the studies of Lauterborn the desired orbit; that is, the application of one of the and Cramer [72] and Lauterborn and Suchla [73] where above targeting techniques is mandatory. complete period-doubling scenarios could be measured If feedback control is possible then the aforemen- in this way of harmonically forced bubbles represented tioned controls work very well; otherwise, the appli- via spectral bifurcation diagrams. Another requirement cation of a non-feedback method (including a stochas- related to the topology is that the transformed trajecto- tic one) is inevitable. The main advantage of the non- ries must be compatible in terms of winding numbers feedback control technique presented in this study is as already discussed in Sect. 8. the ability to target the desired attractor directly, which It must be emphasized that it is possible to combine is not possible with traditional non-feedback controls. the present non-feedback control technique with other For instance, the attractor selection by short pulses [31] methods. For instance, with stochastic control, attrac- (kick the system) introduces randomness to the con- tors with small basins can be annihilated to “simplify” trol; that is, it cannot be guaranteed that the final state the state spaces during the transformation. Moreover, of the system after the kick is the desired attractor. even if feedback control is possible in an application, Moreover, extremely long transients in the presence the method introduced here is still a good candidate for of transient chaos [68] are also possible. With the con- targeting a desired attractor, and after the continuous tinuous transformation of attractors such transients can transformation, the attractor can be stabilized against be avoided. The attractor annihilation techniques by external noise with the feedback control. Similarly, pseudo-periodic forcing [69] or by harmonic pertur- the delay feedback control can stabilize an unstable bation of a state variably or a parameter [32] cannot orbit that replaces the originally stable chaotic attrac- target an attractor directly either. Roughly speaking, tor. After that, the continuous transformation technique both techniques destruct attractors in the order of the can be applied between the resulting coexisting stable size of their basin of attraction with increasing forcing states if necessary and if it is possible. In this way, chaos or perturbation amplitudes. control and targeting an attractor directly is accom- The applicability of the method presented here has plished simultaneously. There are case studies where also requirements/limitations. Firstly, in many cases, the application of two different kinds of controls may there is the requirement that the system should not be be successful when the usage of only one of the control modified appreciably during the control. Therefore, all fails [74]. non-feedback techniques apply only small perturba- Since the control technique is demonstrated numer- tions. It is absolutely not fulfilled in the present method ically only on a specific system (Keller–Miksis equa- as large deviations from the original system parameters tion describing bubble dynamics) and on a very spe- is usually necessary. However, by the end of the control, cific example, many questions arise: Is it possible to the system goes back exactly to the same parameter set “bypass” somehow the winding number compatibility and additional control is no longer required to sustain condition? Is it possible to transform a chaotic attrac- the desired state (attractor). It must be also noted that if tor into a stable periodic orbit and vice versa? Is this the basin of attraction of the targeted attractor is small method applicable to other systems? In order to answer then slow transformation may be necessary. these question correctly, additional detailed analyses Secondly, due to the large parameter variation during are mandatory; however, preliminary conjectures can the control, detailed knowledge on the topology of the be made. periodic orbits in a large parameter space must be avail- Even if a continuous transformation between two able. This is usually a minor problem since the bifur- coexisting attractors does not exist applying a single- 123 Nonlinear oscillators by dual-frequency driving frequency pair, it may be possible to find a transfor- that two attractors which cannot be exchanged by tun- mation by using subsequent frequency pairs one after ing the system parameters may be transformed contin- another, and reaching the desired attractor step-by-step uously into each other by temporarily adding a second hopping onto “intermediate” attractors. The transfor- sinusoidal component to the driving. The process to mation of chaotic attractors into periodic orbits fits into target a desired attractor directly was tested at a rele- the large topic of chaos control [21]. Since chaos con- vant example, the Keller–Miksis equation—a second- trol is possible with the addition of a second sinusoidal order ordinary differential equation describing bubble component to the driving [24,75], the main issue is that dynamics—and visualized via an animation (Online how chaos is annihilated (a smooth inverse Feigenbaum Resource 1). If this technique is proven to be applicable cascade or a sudden crises) during the parameter tuning, to other systems as well, it may become an excellent and what is the resulting periodic orbit which is inter- candidate to accomplish non-feedback control without esting from the winding number compatibility point of long transients and/or randomness. view. It is well known that chaos control is possible Acknowledgements Open access funding provided by Max merely by adjusting the phase shift between the sinu- Planck Society. This paper was supported by the János Bolyai soidal components of the driving [76]; therefore, the Research Scholarship of the Hungarian Academy of Sciences effect of the phase shift must also be thoroughly inves- and the Deutsche Forschungsgemeinschaft (DFG) Grant No. ME tigated. It can also play an important role between the 1645/7-1. transformation of the period-2 and period-3 orbits itself Compliance with ethical standards shown in Fig. 11. The period-2 solutions can always be transformed into one of the period-3 orbits through the Conflict of interest The authors declare that they have no con- blue or the green surfaces. In contrast, if the period-3 flict of interest. attractor is in a phase in time described by the red sur- Open Access This article is distributed under the terms of the face then transformation is not possible in the opposite Creative Commons Attribution 4.0 International License (http:// direction (from period-3 to period-2), see again Fig. 11. creativecommons.org/licenses/by/4.0/), which permits unrestricted It is very likely that the control technique of con- use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, tinuous transformation is possible for other systems provide a link to the Creative Commons license, and indicate if as well. As it was already mentioned, the topology of changes were made. many nonlinear oscillators is composed by the same building blocks. Moreover, in recent publications, a phenomenon called replication of periodic domains is References reported by applying a second sinusoidal component to the driving [40,77], or a single periodic perturbation in 1. Pisarchik, A.N., Feudel, U.: Control of multistability. Phys. case of maps [78,79]. We believe that the replication of Rep. 540(4), 167 (2014) 2. Angeli, D., Ferrell, J.E., Sontag, E.D.: Detection of multista- periodic domains and the branching mechanism intro- bility, bifurcations, and hysteresis in a large class of biolog- duced in Sect. 5.1 have the same dynamical origin. ical positive-feedback systems. Proc. Natl. Acad. Sci. USA Finally, it is worth mentioning that any application 101(7), 1822 (2004) using dual-frequency driving can benefit (at least indi- 3. Ullner, E., Koseska, A., Jürgen, K., Volkov, E., Kantz, H., García-Ojalvo, J.: Multistability of synthetic genetic net- rectly) from the results presented in this study due to works with repressive cell-to-cell communication. Phys. the very detailed investigation. Some possible appli- Rev. E 78, 031904 (2008) cations are acoustic cavitation [80–83], Faraday waves 4. Shiau, Y.H., Peng, Y.F., Hwang, R.R., Hu, C.K.: Multista- [84], stability of traveling beams [85,86]orlaser-driven bility and symmetry breaking in the two-dimensional flow around a square cylinder. Phys. Rev. E 60, 6188 (1999) dissociation of molecules [87]. 5. Hos, ˝ C.J., Champneys, A.R., Paul, K., McNeely, M.: Dynamic behaviour of direct spring loaded pressure relief valves in gas service: II reduced order modelling. J. Loss 10 Conclusion Prevent. Proc. 36, 1 (2015) 6. Ganapathisubramanian, N., Showalter, K.: Bistability, mushrooms, and isolas. J. Chem. Phys. 80(9), 4177 (1984) In the present study, a new non-feedback technique to 7. Yang, L., Dolnik, M., Zhabotinsky, A.M., Epstein, I.R.: control multistability in harmonically driven nonlinear Turing patterns beyond hexagons and stripes. Chaos 16(3), oscillators is introduced. It is based on the observation 037114 (2006) 123 F. Hegedus ˝ et al. 8. Braun, J., Mattia, M.: Attractors and noise: twin drivers of 30. Huisman, J., Weissing, F.: Biodiversity of plankton by decisions and multistability. Neuroimage 52(3), 740 (2010) species oscillations and chaos. Nature 402, 407–410 (1999) 9. Bathiany, S., Claussen, M., Fraedrich, K.: Implications of 31. Chizhevsky, V.N., Grigorieva, E.V., Kashchenko, S.A.: Opti- climate variability for the detection of multiple equilibria mal timing for targeting periodic orbits in a loss-driven CO2 and for rapid transitions in the atmosphere-vegetation sys- laser. Opt. Commun. 133(1), 189 (1997) tem. Clim. Dyn. 38(9), 1775 (2012) 32. Pisarchik, A.N., Goswami, B.K.: Annihilation of one of the 10. Sneppen, K., Mitarai, N.: Multistability with a metastable coexisting attractors in a bistable system. Phys. Rev. Lett. mixed state. Phys. Rev. Lett. 109, 100602 (2012) 84, 1423 (2000) 11. Lauterborn, W., Kurz, T.: Physics of bubble oscillations. 33. Shinbrot, T., Ott, E., Grebogi, C., Yorke, J.A.: Using chaos Rep. Prog. Phys. 73(10), 106501 (2010) to direct trajectories to targets. Phys. Rev. Lett. 65, 3215 12. Lieberman, M.A., Tsang, K.Y.: Transient chaos in dissipa- (1990) tively perturbed, near-integrable hamiltonian systems. Phys. 34. Jiang, Y.: Trajectory selection in multistable systems using Rev. Lett. 55, 908 (1985) periodic drivings. Phys. Lett. A 264(1), 22 (1999) 13. Kolmogorov, A.N.: On the persistence of conditionally peri- 35. Kaneko, K.: Dominance of milnor attractors and noise- odic motions under a small change of the Hamiltonian func- induced selection in a multiattractor system. Phys. Rev. Lett. tion. Dokl. Akad. Nauk. S.S.S.R. 98, 527 (1954). (Russian) 78, 2736 (1997) 14. Feudel, U., Grebogi, C., Hunt, B.R., Yorke, J.A.: Map with 36. Doedel, E.J., Oldeman, B.E., Champneys, A.R., Dercole, more than 100 coexisting low-period periodic attractors. F., Fairgrieve, T.F., Kuznetsov, Y.A., Paffenroth, R., Sand- Phys. Rev. E 54, 71 (1996) stede, B., Wang, X., Zhang, C.: AUTO-07P: Continuation 15. Gavrilov, N.K., Silnikov, L.P.: On three dimensional dynam- and Bifurcation Software for Ordinary Differential Equa- ical systems close to systems with structurally unstable tions. Concordia University, Montreal (2012) homoclinic curve. I. Math. U.S.S.R. Sbornik 17, 467 (1972) 37. Brennen, C.E.: Cavitation and Bubble Dynamics. Oxford 16. Newhouse, S.E.: Diffeomorphisms with infinitely many University Press, New York (1995) sinks. Topology 13(1), 9 (1974) 38. Niemeyer, K.E., Sung, C.J.: Accelerating moderately stiff 17. Wiesenfeld, K., Hadley, P.: Attractor crowding in oscillator chemical kinetics in reactive-flow simulations using GPUs. arrays. Phys. Rev. Lett. 62, 1335 (1989) J. Comput. Phys. 256, 854 (2014) 18. Ikeda, K.: Multiple-valued stationary state and its instability 39. Bonatto, C., Gallas, J.A.C., Ueda, Y.: Chaotic phase similar- of the transmitted light by a ring cavity system. Opt. Com- ities and recurrences in a damped-driven Duffing oscillator. mun. 30(2), 257 (1979) Phys. Rev. E 77(2), 026217 (2008) 19. Pisarchik, A.N.: Dynamical tracking of unstable periodic 40. Medeiros, E.S., de Souza, S.L.T., Medrano-T, R.O., Caldas, orbits. Phys. Lett. A 242(3), 152 (1998) I.L.: Replicate periodic windows in the parameter space of 20. Campos-Mejía, A., Pisarchik, A.N., Arroyo-Almanza, D.A.: driven oscillators. Chaos Solitons Fract. 44(11), 982 (2011) Noise-induced on-off intermittency in mutually coupled 41. de Souza, S.L.T., Lima, A.A., Caldas, I.L., Medrano-T, R.O., semiconductor lasers. Chaos Soliton. Fract. 54, 96 (2013) Guimaraes-Filho, ˝ Z.O.: Self-similarities of periodic struc- 21. Schöll, E., Schuster, H.G.: Handbook of Chaos Control. tures for a discrete model of a two-gene system. Phys. Lett. Wiley, New York (2008) A 376(15), 1290 (2012) 22. Mettin, R., Kurz, T.: Optimized periodic control of chaotic 42. Francke, R.E., Pöschel, T., Gallas, J.A.C.: Zig–zag networks systems. Phys. Lett. A 206(5), 331 (1995) of self-excited periodic oscillations in a tunnel diode and a 23. Mettin, R., Hübler, A., Scheeline, A., Lauterborn, W.: Para- fiber-ring laser. Phys. Rev. E 87(4), 042907 (2013) metric entrainment control of chaotic systems. Phys. Rev. E 43. Medrano-T, R.O., Rocha, R.: The negative side of Chua’s 51, 4065 (1995) circuit parameter space: stability analysis, period-adding, 24. Behnia, S., Sojahrood, A.J., Soltanpoor, W., Jahanbakhsh, basin of attraction metamorphoses, and experimental inves- O.: Suppressing chaotic oscillations of a spherical cavitation tigation. Int. J. Bifurcat. Chaos 24(09), 1430025 (2014) bubble through applying a periodic perturbation. Ultrason. 44. Celestino, A., Manchein, C., Albuquerque, H.A., Beims, Sonochem. 16(4), 502 (2009) M.W.: Stable structures in parameter space and optimal 25. Martínez, P.J., Euzzor, S., Gallas, J.A.C., Meucci, R., ratchet transport. Commun. Nonlinear Sci. Numer. Simul. Chacón, R.: Impulse-Induced Optimum Control of Chaos 19(1), 139 (2014) in Dissipative Driven Systems. ArXiv e-prints (2017) 45. Rech, P.C.: Period-adding structures in the parameter-space 26. Pyragas, K., Lange, F., Letz, T., Parisi, J., Kittel, A.: Stabi- of a driven Josephson junction. Int. J. Bifurc. Chaos 25(12), lization of an unstable steady state in intracavity frequency- 1530035 (2015) doubled lasers. Phys. Rev. E 61, 3721 (2000) 46. da Costa, D.R., Hansen, M., Guarise, G., Medrano-T, R.O., 27. Wieczorek, S., Krauskopf, B., Lenstra, D.: Mechanisms for Leonel, E.D.: The role of extreme orbits in the global organi- multistability in a semiconductor laser with optical injection. zation of periodic regions in parameter space for one dimen- Opt. Commun. 183(1), 215 (2000) sional maps. Phys. Lett. A 380(18), 1610 (2016) 28. Guevara, M., Glass, L., Shrier, A.: Phase locking, period- 47. Freire, J.G., Gallas, M.R., Gallas, J.A.C.: Stability mosaics doubling bifurcations, and irregular dynamics in period- in a forced Brusselator. Eur. Phys. J. Spec. Top. 226(9), 1987 ically stimulated cardiac cells. Science 214(4527), 1350 (2017) (1981) 48. Horstmann, A.C.C., Albuquerque, H.A., Manchein, C.: The 29. Li, E.: Chromatin modification and epigenetic reprogram- effect of temperature on generic stable periodic structures in ming in mammalian development. Nat. Rev. Genet. 3, 662– the parameter space of dissipative relativistic standard map. 673 (2002) Eur. Phys. J. B 90(5), 96 (2017) 123 Nonlinear oscillators by dual-frequency driving 49. Prants, F.G., Rech, P.C.: Complex dynamics of a three- 69. Pecora, L.M., Carroll, T.L.: Pseudoperiodic driving: elim- dimensional continuous-time autonomous system. Math. inating multiple domains of attraction using chaos. Phys. Comput. Simul. 136, 132 (2017) Rev. Lett. 67, 945 (1991) 50. Sewerin, F., Rigopoulos, S.: A methodology for the integra- 70. Scheffczyk, C., Parlitz, U., Kurz, T., Knop, W., Lauterborn, tion of stiff chemical kinetics on GPUs. Combust. Flame W.: Comparison of bifurcation structures of driven dissipa- 162(4), 1375 (2015) tive nonlinear oscillators. Phys. Rev. A 43(12), 6495 (1991) 51. Curtis, N.J., Niemeyer, K.E., Sung, C.J.: An Investigation of 71. Zhang, Y., Zhang, Y., Li, S.: Combination and simultaneous GPU-Based Stiff Chemical Kinetics Integration Methods, resonances of gas bubbles oscillating in liquids under dual- ArXiv e-prints (2016) frequency acoustic excitation. Ultrason. Sonochem. 35, 431 52. Stone, C.P., Alferman, A.T., Niemeyer, K.E.: Accelerating (2017) Finite-Rate Chemical Kinetics with Coprocessors: Com- 72. Lauterborn, W., Cramer, E.: Subharmonic route to chaos paring Vectorization Methods on GPUs, MICs, and CPUs, observed in acoustics. Phys. Rev. Lett. 47(20), 1445 (1981) ArXiv e-prints (2016) 73. Lauterborn, W., Suchla, E.: Bifurcation superstructure in a 53. Hegedus, ˝ F., Hos, ˝ C., Kullmann, L.: Stable period 1, 2 and model of acoustic turbulence. Phys. Rev. Lett. 53(24), 2304 3 structures of the harmonically excited Rayleigh–Plesset (1984) equation applying low ambient pressure. IMA J. Appl. Math. 74. Pisarchik, A.N., Martínez-Zérega, B.E.: Noise-induced 78(6), 1179 (2013) attractor annihilation in the delayed feedback logistic map. 54. Englisch, V., Parlitz, U., Lauterborn, W.: Comparison of Phys. Lett. A 377(42), 3016 (2013) winding-number sequences for symmetric and asymmetric 75. Zhang, Y., Zhang, Y.: Chaotic oscillations of gas bub- oscillatory systems. Phys. Rev. E 92(2), 022907 (2015) bles under dual-frequency acoustic excitation. Ultrason. 55. Fyrillas, M.M., Szeri, A.J.: Dissolution or growth of soluble Sonochem. 40(Part B), 151 (2018) spherical oscillating bubbles. J. Fluid Mech. 277, 381 (1994) 76. Yang, J., Qu, Z., Hu, G.: Duffing equation with two periodic 56. Hos, ˝ C., Champneys, A.R., Kullmann, L.: Bifurcation anal- forcings: the phase effect. Phys. Rev. E 53, 4402 (1996) ysis of surge and rotating stall in the Moore–Greitzer com- 77. Medeiros, E.S., de Souza, S.L.T., Medrano-T, R.O., Caldas, pression system. IMA J. Appl. Math. 68(2), 205 (2003) I.L.: Periodic window arising in the parameter space of an 57. Hegedus, ˝ F., Kullmann, L.: Basins of attraction in a har- impact oscillator. Phys. Lett. A 374(26), 2628 (2010) monically excited spherical bubble model. Period. Polytech. 78. Manchein, C., da Silva, R.M., Beims, M.W.: Proliferation of Mech. Eng. 56(2), 125 (2012) stability in phase and parameter spaces of nonlinear systems. 58. Hos, ˝ C., Champneys, A.R.: Grazing bifurcations and chat- Chaos 27(8), 081101 (2017) ter in a pressure relief valve model. Phys. D 241(22), 2068 79. da Silva, R.M., Manchein, C., Beims, M.W.: Controlling (2012) intermediate dynamics in a family of quadratic maps. Chaos 59. Hegedus, ˝ F.: Stable bubble oscillations beyond Blake’s crit- 27(10), 103101 (2017) ical threshold. Ultrasonics 54(4), 1113 (2014) 80. Holzfuss, J., Rüggeberg, M., Mettin, R.: Boosting sonolu- 60. Hegedus, ˝ F., Klapcsik, K.: The effect of high viscosity on minescence. Phys. Rev. Lett. 81, 1961 (1998) the collapse-like chaotic and regular periodic oscillations of 81. Krefting, D., Mettin, R., Lauterborn, W.: Two-frequency a harmonically excited gas bubble. Ultrason. Sonochem. 27, driven single-bubble sonoluminescence. J. Acoust. Soc. Am. 153 (2015) 112(5), 1918 (2002) 61. Hegedus, ˝ F.: Topological analysis of the periodic structures 82. Yasuda, K., Torii, T., Yasui, K., Iida, Y., Tuziuti, T., Naka- in a harmonically driven bubble oscillator near Blake’s criti- mura, M., Asakura, Y.: Enhancement of sonochemical reac- cal threshold: infinite sequence of two-sided Farey ordering tion of terephthalate ion by superposition of ultrasonic fields trees. Phys. Lett. A 380(9–10), 1012 (2016) of various frequencies. Ultrason. Sonochem. 14(6), 699 62. Parlitz, U., Lauterborn, W.: Resonances and torsion numbers (2007) of driven dissipative nonlinear oscillators. Z. Naturforsch. A 83. Merouani, S., Hamdaoui, O., Rezgui, Y., Guemini, M.: Sen- 41(4), 605 (1986) sitivity of free radicals production in acoustically driven bub- 63. Parlitz, U., Lauterborn, W.: Period-doubling cascades and ble to the ultrasonic frequency and nature of dissolved gases. devil’s staircases of the driven van der Pol oscillator. Phys. Ultrason. Sonochem. 22, 41 (2015) Rev. A 36(3), 1428 (1987) 84. Batson, W., Zoueshtiagh, F., Narayanan, R.: Two-frequency 64. Lai, Y.C.: Driving trajectories to a desirable attractor by excitation of single-mode Faraday waves. J. Fluid Mech. using small control. Phys. Lett. A 221(6), 375 (1996) 764, 538 (2015) 65. Gadaleta, S., Dangelmayr, G.: Optimal chaos control 85. Sahoo, B., Panda, L.N., Pohit, G.: Two-frequency parametric through reinforcement learning. Chaos 9(3), 775 (1999) excitation and internal resonance of a moving viscoelastic 66. Martínez-Zérega, B.E., Pisarchik, A.N., Tsimring, L.S.: beam. Nonlinear Dyn. 82(4), 1721 (2015) Using periodic modulation to control coexisting attractors 86. Sahoo, B., Panda, L.N., Pohit, G.: Combination, principal induced by delayed feedback. Phys. Lett. A 318(1–2), 102 parametric and internal resonances of an accelerating beam (2003) under two frequency parametric excitation. Int. J. Nonlinear 67. Poon, L., Grebogi, C.: Controlling complexity. Phys. Rev. Mech. 78(Supplement C), 35 (2016) Lett. 75, 4023 (1995) 87. Huang, S., Chandre, C., Uzer, T.: Bifurcations as dis- 68. Lai, Y.C., Tél, T.: Transient Chaos. Springer, New York sociation mechanism in bichromatically driven diatomic (2010) molecules. J. Chem. Phys. 128(17), 174105 (2008)

Journal

Nonlinear DynamicsSpringer Journals

Published: May 29, 2018

References

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off