Multi-stability in Doubochinski’s Pendulum

Multi-stability in Doubochinski’s Pendulum

Yao Luo School of Physics, Nanjing University, Nanjing, P.R China, 210093    Wenkai Fan Kuang Yaming Honors School, Nanjing University, Nanjing, P.R China, 210093 Department of Physics, Duke University    Chenghao Feng Kuang Yaming Honors School, Nanjing University, Nanjing, P.R China, 210093 Department of Electrical and Computer Engineering, The University of Texas at Austin.    Sihui Wang School of Physics, Nanjing University, Nanjing, P.R China, 210093    Yinlong Wang School of Physics, Nanjing University, Nanjing, P.R China, 210093
July 26, 2019

The widespread phenomena of multistability is a problem involving rich dynamics to be explored. In this paper, we study the multistability of a generalized nonlinear forcing oscillator excited by . We take Doubochinski’s Pendulum as an example. The so-called ”amplitude quantization”, i.e., the multiple discrete periodical solutions, is identified as self-adaptive subharmonic resonance in response to nonlinear feeding. The subharmonic resonance frequency is found related to the symmetry of the driving force: odd subharmonic resonance occurs under even symmetric driving force and vice versa. We solve the multiple periodical solutions and investigate the transition and competition between these multi-stable modes via frequency response curves and Poincare maps. We find irreversible transition between the multistable modes and propose a multistability control strategy.

Valid PACS appear here
preprint: APS/123-QED

I Introduction

The phenomenon of multistability (coexistence of different attractors for a given set of parameters) is widespread in physics maurer1980effect (); brun1985observation (); gibbs2012optical (), chemistry aguda1987bistability (); wilhelm2009smallest (), biology angeli2004detection (); ozbudak2004multistability (); freyer2011biophysical () and climate systems robinson2012multistability (); freire2008multistability () and many other fields of science and nature. One important feature of such systems of high order of complexity are their extreme sensitivity to perturbation and initial conditionspisarchik2014control (). Generalized multistability, firstly distinguished from ordinary bistability extensively exist in nonlinear systems, was found in laser physicsarecchi1982experimental ().

The mechanisms of multistability include coupling maistrenko2007multistability (), delayed feedback balanov2005delayed (); foss1996multistability (),parametrical forcing varma1993quadratic () and dissipationlieberman1985transient (). These mechanisms usually do not act independently in a system. In many cases, the feedback and dissipation often coexist. For example, a mechanical rotor with weak damping and self-adaptive driving term possess over one hundred attractorsfeudel1996map ().

Compared to bistability phenomenon, the study of generalized multistability is far from being complete. The phenomena and mechanisms of multistability involve rich dynamics to be explored. Here we study the multiple stability behaviors in a classical mechanical system, the Doubochinski’s pendulumtennenbaum2006amplitude (). We reveal that it’s mysterious ”amplitude quantization”, initially studied as ”macroscopic quantum behavior”tennenbaum2006amplitude () is, in fact, the self-adaptive behavior in response to nonlinear feeding.   We consider a generalized model of an oscillator excited by periodic force where the feeding function is a nonlinear function of displacement. This model containing restoring force, dissipation and nonlinear feeding can be regarded as generalization of parametric oscillation. For a linear parametric oscillator governed by Mathieu equation, subharmonic resonance occurs when the driving frequency is near twice the natural frequency of an oscillator nayfeh2008nonlinear (). To understand the mechanism of multistability, we express the feeding function as series of polynomial and solved the amplitudes and frequency respond curves analytically. The origin of multiple solutions and strong self-adaptivity can be explained. Moreover, the even-odd correspondence between the symmetry of the driving force and subharmonic resonance frequencies are presented numerically. To investigate the transition and competition between the multi-stable modes, we present the frequency response curves and Poincare maps near and . The transition between multistable states is more complex and intriguing than in bistability. We find irreversible transition between the multistable modes and propose a multistability control strategy. In contrast, the frequency response for bistability states forms a closed hysteresis loop rather than an open route.


Ii Theoretical Model

ii.1 Dynamic Equation

A general form of nonlinearly excited oscillator can be described by the following dynamic equation


where is the displacement, is the linear damping coefficient, is the restoring force of the system, is a nonlinear periodic driving force and is the driving frequency. is the displacement-dependent amplitude of the driving force, called feeding function here. There are various physical systems described by equations of similar forms. For example, when contains a positive linear term and a cubic term, and the feeding function is a constant, Eq.1 becomes Duffing’s equation kovacic2011duffing (). A nonlinear parametric oscillator is described by an equation similar to Eq.1, where is proportional to displacement. rhoads2006generalized ().

When applied to a pendulum the restoring force can be written as


where is the pendulum’s natural frequency and is the angular displacement.

The actual form of the feeding function is cumbersome, so we approximate the feeding function in terms of polynomial of to do analytical calculation. Specifically, if is an even symmetric function of , the polynomial approximation only contains even-degree terms,


If is an odd symmetric function of , it only contains odd-degree terms.


In many cases, the feeding function is considerable in an active zone and becomes negligible elsewhere, such a pendulum is called a ”kick-excited” one. In Ref. damgov2000discrete (), a symmetric shaped feeding function is considered. Here, the polynomial series is truncated after the lowest few terms to model the feeding function. The coefficients are determined by fitting of a particular physical model using above polynomials. Since the polynomial usually does not converge to the practical feeding function for arbitrary displacement, the polynomial approximation only applies to the range of oscillation or the ”active zone” of the feeding function unnecessarily being small.

Figure 1: a,b, Schematic diagrams of the Doubochinski’s pendulum. The electromagnet is vertically(a) or horizontally(b) placed beneath the pendulum. c,d, Feeding function of Doubochinski’s pendulum for vertically(c) or horizontally(d) placed electromagnet.

ii.2 Feeding function for a magnetic pendulum

Later on, we consider a practical model of a nonlinearly excited oscillator, the Doubochinski’s pendulum. The pendulum consists of a light rigid pendulum with a small magnet at the free end. An AC powered electromagnet is vertically or horizontally placed beneath the pendulum, see Figs.1 and 1. The magnet pole’s direction should be adjusted along the local geomagnetic field to eliminate the torque. The dynamic equation for the pendulum is


where is the moment of inertia of the pendulum, is the torque applied to the magnet. We can write the dynamical equation in a dimensionless form.


where , and are defined as


The driving term arises from the alternating magnetic field. The torque can be evaluated by treating the magnet as a magnetic dipole, thus,


where is the magnetic induction of the field, and is the magnetic moment of the magnet. Here we calculate the feeding function in two typical configurations, parameters for calculation are given in table 1, involving natural frequency of the pendulum , moment of inertia , length of the pendulum , radius of the coils , number of turns , height of the electromagnet coil , distance between the coil and the magnet , magnetic momentum of the magnet and the current . These physical parameters will be used for numerical calculation. The coefficients of the polynomial , , in Eq.3 fitting to the feeding function given by above parameters are also listed, which will be used in analytical calculation. As shown in Figs.1 and 1, when the electromagnetic coil is placed vertically, has even symmetry; when the coil is placed horizontally, has odd symmetry.

Parameters :
5.13 0.01 0.456 0.042 220
0.02 0.021 1.36 0.5 0.0156
fitting coefficients of the polynomial :
0.055 -11.6 300
Table 1: The parameters used to calculate the feeding function and the fitting coefficients of the polynomial.

Iii Dynamic Behaviors

iii.1 Discrete Periodic Orbits

The pendulum’s motion is governed by Eq.6 and Eq.9. We solve the equations when the electromagnet is vertically placed, using the feeding function as shown in Fig.1 and Runge-Kutta 4th order method in c++.

In Fig.2, we find that there are several stationary oscillation modes near triple () and quintuple of the natural frequency ().Figs.2(a) and 2(b) are the time evolutions of typical periodic orbits. The final oscillation modes are dependent on initial conditions. Then we sample the initial condition in phase space of , in search of stable oscillation orbits. All the trajectories converge to several limit cycles, as shown in Figs.2(c) and 2(d). The small cycles in the middle are the linear solutions where the pendulum oscillates at the frequencies of the driving force. The larger cycles are nonlinear solutions where the pendulum oscillates near the natural frequency. The insets of Figs.2(c) and 2(d) show the dual stationary solutions corresponding to the larger limit cycles. Though the two limit cycles are adjacent in phase diagram, they are different in phase. From Fig.2, we also see that for the large stationary solutions, the pendulum approximately performs harmonic oscillation.

These results are consistent with experimental phenomenon njubook () that pendulum released at small angles only does small oscillations in the ”linear region”; and for larger releasing angles, the amplitude varies slowly and ends up with stable harmonic oscillations on the large amplitude orbits or the small linear orbit.

Figure 2: a,b, Multiple stationary solutions of Doubochinski’s pendulum for driving frequency . The final oscillation modes vary with initial conditions. c,d, Limit cycles for respectively. The small cycles in the middle are the linear solutions. The larger cycles which appear in pairs are nonlinear solutions where the pendulum oscillates near the natural frequency. The insets show stationary solutions of the corresponding limit cycles in the same color.

To find the multiple periodic solutions analytically, we use a polynomial approximation of the feeding function Eq.3 and truncate it after the 4th order.

In Eq.3, when , the leading term is a constant . The pendulum only does regular linear forced oscillation at the driving frequency. At larger angles, the effect of nonlinearity becomes significant. When the pendulum performs large amplitude oscillation under driving frequency , the zeroth order solution of Eq.6 should be written as


In order to capture the frequency response of the system, a factor is added to the triple natural frequency


Substituting Eq.10 into Eq.6 and averaging the amplitude and phase change over one period, we get the time derivatives of amplitude and phase,


For simplicity, we set the damping coefficient . These two equations approximately describe the evolution trajectories in space. The stationary solutions are the fixed points in space, which satisfy


The solutions are


in which .

Figure 3: The analytical result of amplitude frequency response curves for . The two black solid curves denote the larger stable roots, the blue dashed the unstable root. The linear forced oscillation solution is the solid line near the bottom.

The term can be ignored in further analysis because it is much smaller than and , see table.1. When , the amplitudes solved from Eq15 and 16 are the same,


Eq.17 indicates that nontrivial stationary solutions exist when and are different in sign. Here multi-stability is achieved without dissipation when the polynomial of the feeding function contains more than one nonlinear terms so that energy balance can be realized. Feeding function truncated after the 4th order can produce resonance at . For higher subharmonic frequencies, feeding function with higher order terms should be taken into consideration.

The oscillation amplitudes in Eq.17 are independent of the amplitude of the driving force, since the , are proportional to the amplitude of the driving force and the ratio will not change. Such strong self-adaptivity is verified in Ref.damgov2000discrete (); njubook (), in which the amplitudes almost do not change with the driving force within a very large range. Actually, as amplitude of the driving force exceeds a critical value, the oscillator will undergo a series of symmetry breaking oscillation modes with multiple period and finally becomes chaotic. The period-3 bifurcations in a system with shaped feeding function reported by Domgov damgov2000discrete () is a similar example.

The frequency response curves shown in Fig.3 can be solved from Eq.16 when . The nonlinear solutions contain one smaller unstable root, denote by blue dotted curve, and a pair of larger stable roots, denote by black solid curves. The linear forced oscillation solution is also plotted, see the solid line near the bottom in Fig.3.

Figure 4: The max-amplitude frequency response diagram. The red (blue) curves are for even (odd) feeding function. Subharmonic resonance occurs for in response to even (odd) feeding function. The frequency response curve for linear forced oscillation is not plotted.

iii.2 Symmetry of the Driving Force and Frequency Response

To clarify how symmetry of the feeding function influences the subharmonic frequencies, we calculate the overall frequency response diagram for even symmetric and odd symmetric . In calculation, we consider damping and using Eq.9 and parameters listed in table.1 to model the feeding function, the frequency response diagram can be solved numerically. Since multiple periodic oscillation modes coexist when subharmonic resonance occurs, we only show maximum amplitudes on the frequency response diagram.

In Fig.4, red curves correspond to even symmetric feeding function, and blue curves for odd symmetric feeding function. We find that subharmonic resonance occurs for () in response to even (odd) symmetric feeding function. Since an asymmetric function can be decomposed into an odd function and even function, we can infer that if the electromagnet is placed inclined, subharmonic resonance could occur for .

Fig.5(e) shows the frequency response diagram when . Similar to analytical result in Fig.3, it has a small linear solution, denoted as mode , and a pair of intersecting branches denoted as mode and mode for . and stand for right and left; , and denote different branches. Here the unstable solutions are not presented because they are unavailable via numerical integral of Eq.6.

Figure 5: a,b,c,d, Evolution of attraction basins in Poincare maps with varying . The black dots are FPs for each basin. Attraction basin for each FP is denoted with a color. e, Frequency response curves for . The arrows beside the frequency response curves indicate how oscillation modes evolve with driving frequency.

Fig.6(e) is the frequency response diagram for . The five stable orbits in Fig.2(d) correspond to the one linear branch named and two pairs of overlapped branches, named as .

iii.3 Modes Competition

A Duffing oscillator with cubic nonlinearity in restoring forcekovacic2011duffing () exhibits jump phenomena and hysteresis between bistable states due to hardening or softening resonance. Meanwhile, the frequency response curves of a nonlinear parametric oscillator with bistability exhibit mixed feature of softening and hardeningrhoads2006generalized (). Here, more complex jump phenomenon, induced by modes competition among multiple attractors, is studied.

In order to study the transition behavior among the attractors, the frequency response graph is calculated when excitation frequency increases or decreases quasi-statically. Transition happens when one orbit loses stability and the oscillator jumps to another orbit, as the control parameter, i.e. the driving frequency, changing quasi-statically. In Fig.5(e), the arrows beside the frequency response curves indicate how transition happens. The frequency condition when all stationary modes , and coexist is . When , modes and are stable, and when , modes and are stable. As frequency increases, the oscillator in mode jumps to mode at . Unexpectedly, as frequency decreases, transition from mode to mode happens at .

Such irreversible transition behavior can be described from the evolution of attraction basins in Poincare maps. Figs.5(a) to 5(d) are the attraction basins when respectively. In each figure, the small black dots are fixed points (FPs). Each FP is surrounded by its basin of attraction indicated with color. The red, green and purple areas denote the attraction basins for modes , and respectively. Mode (the trivial solution) leaves one fixed point near the center in Figs.5(a) to 5(d). Meanwhile, each large-amplitude mode has 3 fixed points because the oscillation periods of modes and are triple of excitation period.

In Figs.5(b), 5(c) and 5(d), as the excitation frequency increases from to 3.00 and 3.040, the attraction basin of shrinks quickly and vanished below 3.040. In Fig.5(c) just below the transition frequency, the vanishing attraction basins of surrounded by those of will merged into those of and transition from to happens.

In Figs.5(b) and 5(a), in vicinity of the critical frequency , the basins of are surrounded by those of mode (the red area). As the excitation frequency decreases from 2.945 to 2.940, the basins of mode are replaced by those of mode so that transition from mode to happens.

Figure 6: a,b,c,d, Evolution of the attraction basins in Poincare maps for . The black dots are FPs for each basin. Attraction basin for each FP is denoted with a color.e, Frequency response curves for . The arrows beside the frequency response curves indicate how oscillation modes evolve with driving frequency.

Additionally, the frequency response diagram and Poincare maps for are presented, see Fig.6. is the lower critical frequency for mode , and is the upper critical frequency for mode . In Figs.6(a) to 6(d), the red, orange, green, purple and blue areas are the attraction basins of modes , , , and respectively. We mainly focus on a typical transition starting from model and mode below.

Fig.6(e) shows that transition among different modes is also irreversible. When increasing the driving frequency, oscillator at mode will jump to mode at . In Fig.6(c), the attraction basins of mode are surrounded by those of below , and the basins of will merge into the basin as driving frequency exceeds , see Fig.6(d).

When decreasing the driving frequency, oscillator at mode jumps to mode at , rather than . In Fig.6(b) and Fig.6(a), the attraction basins of mode will merge into the basin of mode below .

From the transition sequence described in Fig.5(e) and Fig.6(e), we propose a control strategy among the oscillation modes by adiabatically changing the driving frequency. For example, in Fig.5(e) starting from mode , the oscillator can be induced to mode by increasing the driving frequency, then to mode by decreasing driving frequency. However, the transition from mode to higher oscillation modes are forbidden. It means that an oscillator initially performs linear oscillation in mode , cannot be driven to higher orbits merely by changing the driving frequency.

Iv Conclusion

The phenomena and mechanisms of multistability involve rich dynamics to be explored. In this paper, we investigate the multistability of a generalized nonlinear forcing oscillator excited by . We take Doubochinski’s Pendulum as an example. The multistability mechanism in our system is nonlinear parametric resonance. We express the energy feeding function into polynomial and find multi-stability can be achieved without dissipation when the polynomial contains more than one nonlinear terms so that energy balance can be realized. We also find the subharmonic resonance frequency conditions and demonstrate the even-odd correspondence between the symmetry of the feeding function and subharmonic resonance frequencies.

We present the frequency response diagrams and Poincare maps near subharmonic frequencies and . The frequency response diagram contains one linear response branch and nonlinear branches appearing in pairs. We find irreversible transition phenomenon between the multistable modes. Control of the multi-stability can be realized by changing the excitation frequency adiabatically.


  • (1) J. Maurer and A. Libchaber, Effect of the prandtl number on the onset of turbulence in liquid 4he, Journal de Physique lettres. 41(21), 515–518 (1980).
  • (2) E. Brun, B. Derighetti, D. Meier, R. Holzner, and M. Ravani, Observation of order and chaos in a nuclear spin–flip laser, JOSA B. 2(1), 156–167 (1985).
  • (3) H. Gibbs, Optical bistability: controlling light with light. Elsevier (2012).
  • (4) B. D. Aguda and B. L. Clarke, Bistability in chemical reaction networks: theory and application to the peroxidase–oxidase reaction, The Journal of chemical physics. 87(6), 3461–3470 (1987).
  • (5) T. Wilhelm, The smallest chemical reaction system with bistability, BMC systems biology. 3(1), 90 (2009).
  • (6) D. Angeli, J. E. Ferrell, and E. D. Sontag, Detection of multistability, bifurcations, and hysteresis in a large class of biological positive-feedback systems, Proceedings of the National Academy of Sciences. 101(7), 1822–1827 (2004).
  • (7) E. M. Ozbudak, M. Thattai, H. N. Lim, B. I. Shraiman, and A. Van Oudenaarden, Multistability in the lactose utilization network of escherichia coli, Nature. 427(6976), 737 (2004).
  • (8) F. Freyer, J. A. Roberts, R. Becker, P. A. Robinson, P. Ritter, and M. Breakspear, Biophysical mechanisms of multistability in resting-state cortical rhythms, Journal of Neuroscience. 31(17), 6353–6361 (2011).
  • (9) A. Robinson, R. Calov, and A. Ganopolski, Multistability and critical thresholds of the greenland ice sheet, Nature Climate Change. 2(6), 429 (2012).
  • (10) J. G. Freire, C. Bonatto, C. C. DaCamara, and J. A. Gallas, Multistability, phase diagrams, and intransitivity in the lorenz-84 low-order atmospheric circulation model, Chaos: An Interdisciplinary Journal of Nonlinear Science. 18(3), 033121 (2008).
  • (11) A. N. Pisarchik and U. Feudel, Control of multistability, Physics Reports. 540(4), 167–218 (2014).
  • (12) F. Arecchi, R. Meucci, G. Puccioni, and J. Tredicce, Experimental evidence of subharmonic bifurcations, multistability, and turbulence in a q-switched gas laser, Physical Review Letters. 49(17), 1217 (1982).
  • (13) Y. L. Maistrenko, B. Lysyansky, C. Hauptmann, O. Burylko, and P. A. Tass, Multistability in the kuramoto model with synaptic plasticity, Physical Review E. 75(6), 066207 (2007).
  • (14) A. G. Balanov, N. B. Janson, and E. Schöll, Delayed feedback control of chaos: Bifurcation analysis, Physical Review E. 71(1), 016222 (2005).
  • (15) J. Foss, A. Longtin, B. Mensour, and J. Milton, Multistability and delayed recurrent loops, Physical Review Letters. 76(4), 708 (1996).
  • (16) V. Varma et al., Quadratic map modulated by additive periodic forcing, Physical Review E. 48(3), 1670 (1993).
  • (17) M. A. Lieberman and K. Y. Tsang, Transient chaos in dissipatively perturbed, near-integrable hamiltonian systems, Physical review letters. 55(9), 908 (1985).
  • (18) U. Feudel, C. Grebogi, B. R. Hunt, and J. A. Yorke, Map with more than 100 coexisting low-period periodic attractors, Physical Review E. 54(1), 71 (1996).
  • (19) J. Tennenbaum, Amplitude quantization as an elementary property of macroscopic vibrating systems, 21ST CENTURY SCIENCE AND TECHNOLOGY. 18(4), 50 (2006).
  • (20) A. H. Nayfeh and D. T. Mook, Nonlinear oscillations. John Wiley & Sons (2008).
  • (21) I. Kovacic and M. J. Brennan, The Duffing equation: nonlinear oscillators and their behaviour. John Wiley & Sons (2011).
  • (22) J. F. Rhoads, S. W. Shaw, K. L. Turner, J. Moehlis, B. E. DeMartini, and W. Zhang, Generalized parametric resonance in electrostatically actuated microelectromechanical oscillators, Journal of Sound and Vibration. 296(4-5), 797–829 (2006).
  • (23) V. Damgov and I. Popov, “discrete” oscillations and multiple attractors in kick-excited systems, Discrete Dynamics in Nature and Society. 4(2), 99–124 (2000).
  • (24) W. G. Sihui Wang, International Young Physicists’ Tournament: Problems and Solutions 2015. WorldScientific (2018).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description