Simple analytical expressions are suggested for transition curves that separate, in the Ince–Strutt diagram, different types of solutions to the famous Mathieu equation. The derivations of these expressions in this paper rely on physically meaningful periodic solutions describing various regular motions of a familiar nonlinear mechanical system—a rigid planar pendulum with a vertically oscillating pivot. The paper is accompanied by a relevant simulation program.

## I. MATHIEU EQUATION IN PHYSICS AND ENGINEERING

The famous Mathieu equation is an ordinary second-order linear homogeneous differential equation with periodic coefficients that was first introduced by French mathematician E. Leonard Mathieu, in his “Memoir on the vibrations of an elliptic membrane” in 1868.^{1} In the standard notation this equation is usually written as

where the real constants *a* and *q* are often referred to as the characteristic number and parameter, respectively. The Mathieu equation and its solutions^{2} are too well known to require a detailed introduction. This equation is encountered in many different issues in physics, engineering and industry, including the stability of floating ships and railroad trains, the motion of charged particles in electromagnetic Paul traps, the theory of resonant inertial sensors, and many other problems (see, for example, Ref. 3). Its solutions govern the behavior of physical systems of great diversity, and have accordingly been the subjects of a vast number of investigations. There exists a comprehensive literature about the heavily studied Mathieu equation. Reference 4 contains an extensive bibliography of earlier papers and books on the subject.

We shall write Mathieu equation using a slightly different notation from that in Eq. (1), a notation which will be more convenient for our further discussion, namely, as

For the great majority of problems associated with the Mathieu equation, the crucial issue is the determination of conditions in which solutions to this equation remain bounded in the course of time, or grow indefinitely. The answer to this question is given by the well-known Ince–Strutt diagram (Fig. 1), which shows the transition curves in the parameters plane. These curves divide the plane $(k,\u2009m)$ into regions that correspond to unbounded (unstable) solutions (shaded areas in Fig. 1) and stable motions. Stability charts of the Mathieu equation have been investigated using a variety of mathematical methods and concepts such as Lyapunov exponents, Poincaré maps, Floquet transformations, and various perturbation techniques.

There are hundreds of papers and book chapters, in which various algorithms and numerical schemes are suggested for calculating the stability diagram for the Mathieu equation. In this paper, we use an alternative physically transparent approach, which allowed us to find, for the first time, fairly simple analytical formulas for expressing the instability boundaries, in contrast to the numerical algorithms described in the literature.

We note that our formulas give the stability boundaries in a finite form, while all other known analytical derivations such as the two-variable expansion method (see, for example, Ref. 5) deliver an expression for a particular transition curve in the form of a power series expansion.

Our derivation of the stability diagram depends to a large extent on the use of the physical meaning of periodic motions in a familiar mechanical system—the rigid planar pendulum with a vertically oscillating pivot, which is described by a nonlinear generalization of the Mathieu equation. We analyze the relevant periodic solutions to this equation in terms of physics. Numerical simulations show that such solutions are characterized by rather simple spectral composition. This allows us, using a transition to the limiting case of the corresponding linear system described by the Mathieu equation, to express the transition curves of the Ince–Strutt diagram (Fig. 1) in a finite form by the following simple expressions (which are approximate, though sufficiently accurate for applications):

- Curve 2—the upper boundary of dynamic stabilization of the pendulum in the inverted position (at
*k*< 0), and the left (sine-like) boundary of the principal parametric resonance (at $0<k<1/4$), see Sec. III:(4)$m2(k)=14((9\u22124k)(13\u221220k)\u2212(9\u22124k)),\u2003(k<1/4).$ - Curve 3—the right (cosine-like) boundary of the principal parametric resonance (emerging from
*m*= 0 at subharmonic frequency $\omega /2$ of the drive frequency*ω*), see Sec. III:(5)$m3(k)=14(9\u22124k\u2213(9\u22124k)(13\u221220k)),\u2003(1/4<k<13/20).$ - Curve 4—the left (sine-like) boundary of the second-order parametric resonance (at the drive frequency
*ω*), see Sec. IV:(6)$m4(k)=2(k\u22121)(k\u22124)(k\u22129)k\u22125,\u2003(13/20<k<1).$ - Curve 5—the right (cosine-like) boundary of the second-order parametric resonance, see Sec. IV:(7)$m5(k)=2k(k\u22121)(k\u22124)3k\u22128,\u2003(k>1).$

The transition curves of the Ince–Strutt stability diagram in Fig. 1 are plotted with the help of analytical formulas (3)–(7), whose derivation is based on relevant periodic solutions to the Mathieu equation. The physical system that is used for interpretation of these solutions is described in Sec. II of this paper; the derivation of expressions (3)–(7) is discussed in Secs. III–V.

## II. THE PHYSICAL SYSTEM

Mathieu equation governs the response of many physical systems to a sinusoidal parametric forcing. The importance of parametric excitation for numerous applications has led to a wealth of literature on the theoretical analysis and experimental investigation of parametrically excited systems (see, for example, Ref. 6 and references therein). In a great majority of papers and monographs, the phenomena associated with parametric excitation and instability are explained in terms of the theory of differential equations with periodic coefficients (Floquet theory, infinite Hill determinants, continued fractions, etc.; see, for example, Refs. 7 and 8). Most of such papers are predominantly mathematical in nature and actually give little insight into the investigated phenomena whose physical sense, as a rule, is buried deeply in hard mathematics. It is not always easy to understand the physics of the parametric excitation from the general mathematical theory, which could turn out to be abstract and very complicated for physicists and engineers. In this paper, on the contrary, we present a simple treatment of the transition curves of the Ince–Strutt diagram appealing to their physical meaning in order to make the corresponding mathematics more transparent.

To find the transition curves of the Ince–Strutt diagram, we rely on the physically understandable familiar behavior of a well-known simple physical system. Namely, we refer to a rigid planar pendulum with the vertically oscillating pivot, which is described by a nonlinear generalization of the Mathieu equation.

A fascinating feature in the behavior of a simple rigid pendulum whose suspension point is forced to vibrate at a high frequency along the vertical line is the dynamic stabilization of the inverted position. When the frequency and the amplitude of these vibrations are large enough, the inverted pendulum shows no tendency to turn down. Moreover, at moderate deviations from the inverted vertical position, the pendulum tends to return to it. When deviated, the pendulum executes relatively slow oscillations about the vertical line on the background of rapid oscillations of the suspension point. Such extraordinary behavior of the pendulum was physically explained and investigated experimentally in detail by Pjotr Kapitza,^{9} and the corresponding physical device is now widely known as “Kapitza's pendulum.” Conditions of dynamic stabilization of the inverted position give equations for curves 1 and 2 of the Ince–Strutt diagram (Fig. 1). Curves 2 and 3 are bounding the “tongue” of the principal parametric resonance in an ordinary (hanging down) pendulum with vertically oscillating pivot, while curves 4 and 5 correspond to the parametric resonance of the second order.

We consider for simplicity a light rigid rod of length *l* with a heavy small bob of mass *M* on its end and assume that all the mass of the pendulum is concentrated in the bob. The force of gravity $Mg$ (here $g$ is the free-fall acceleration) provides a restoring torque $\u2212Mgl\u2009sin\u2009\phi $ whose value is proportional to the sine of angular deflection $\phi $ of the pendulum from the lower equilibrium position. With the suspension point at rest, this torque makes the deviated pendulum swing about the lower stable equilibrium position. When the axis of the pendulum is forced to move with acceleration along the vertical line, it is convenient to analyze the motion in the non-inertial reference frame associated with this axis. Due to the acceleration of this frame of reference, the pseudo force of inertia $\u2212Mz\xa8$ is exerted on the pendulum, where *z*(*t*) is the momentary vertical coordinate of the axis. The torque $\u2212Mz\xa8l\u2009sin\u2009\phi (t)$ of this time-dependent force must be added to the torque of the gravitational force.

We assume that the axis is forced to execute a given harmonic oscillation along the vertical line at frequency *ω* and amplitude *a*: $z(t)=a\u2009cos\u2009\omega t.$ The force of inertia $Fin(t)$, exerted on the bob in the non-inertial frame of reference, also has the same sinusoidal dependence on time

This force of inertia is directed upward during the time intervals for which *z*(*t*) > 0, i.e., when the axis is located above the middle point of its oscillations. Therefore, during the corresponding half-period of the oscillation of the pivot, this additional force is equivalent to some weakening of the force of gravity. During the other half-period, the axis is below its middle position (*z*(*t*) < 0), and the action of this additional force is equivalent to some strengthening of the gravitational force.

The graphs of the time history and the phase trajectories for different modes of the pendulum, presented later in this paper to illustrate the transition curves of the Ince–Strutt diagram, are obtained with the help of the simulation program “Pendulum with the vertically driven pivot,” which is a part of the software package^{10} “Nonlinear Oscillations” developed by the author. The simulations are based on a numerical integration of the exact differential equation for the angular deflection $\phi (t)$ of the pendulum with the vertically oscillating suspension point. This equation includes, besides the torque $\u2212Mgl\u2009sin\u2009\phi $ of the gravitational force *Mg*, the instantaneous torque of the force of inertia $Fin(t)$, which depends explicitly on time *t*,

The second term in Eq. (9) represents the frictional torque, which is assumed in this model to be proportional to the momentary value of the angular velocity $\phi \u0307$. The damping constant *γ* in this term is inversely proportional to quality factor *Q*, which is conventionally used to characterize damping of small natural oscillations under viscous friction: $Q=\omega 0/2\gamma $, where $\omega 0=g/l$ is the frequency of infinitely small natural oscillations of the pendulum in the absence of pivot oscillations. For small $\phi $ values ($\phi \u226a1$), and in the absence of friction, Eq. (9) reduces to the Mathieu equation, Eq. (2), in which *m* is the normalized amplitude of the pivot ($m=a/l$), and parameter *k* is defined by the following expression:

This dimensionless parameter *k* (inverse normalized drive frequency squared), being physically less meaningful than $\omega /\omega 0$, is nevertheless more convenient for further investigation, because the curves in the Ince–Strutt diagram are usually expressed in terms of *k*. Negative values of *k* correspond to negative *g* values, that is, to upward “gravity.” Negative *k*-values can be formally assigned to the inverted pendulum in the ordinary (directed downward) gravitational field.

## III. THE PRINCIPAL PARAMETRIC RESONANCE: CURVES 2–3 OF THE INCE–STRUTT DIAGRAM

The principal parametric resonance in a pendulum with a vertically oscillating pivot occurs (at small driving amplitudes $m=a/l\u226a1$) if the drive frequency *ω* is approximately twice the natural frequency *ω*_{0} of the pendulum: $\omega \u22482\omega 0$. At a finite amplitude *a* of the pivot, parametric excitation occurs in an interval of frequencies near $2\omega 0$. The greater the value of $m=a/l$, the wider this interval is. In an idealized linear system, within this interval of parametric instability the amplitude of initially small oscillations grows without limit even in the presence of friction, if *m* exceeds some threshold value.

In a real system like the pendulum, the growth of the amplitude at parametric resonance is restricted by nonlinear effects. Such behavior is illustrated in Fig. 2. As the amplitude grows, the natural period of the pendulum becomes longer. If conditions for parametric excitation are fulfilled at small oscillations, and the amplitude is growing, the conditions of resonance become violated at large amplitudes. The drive drifts out of phase with the pendulum, and the phase relations change gradually to those favorable for gradual reduction of the amplitude. The natural period becomes shorter, and conditions for the growth of the amplitude restore. Oscillations acquire the character of beats. Due to friction, these transient beats gradually fade, and a steady-state regular period-2 regime (limit cycle) is established. Such a period-2 regime, characterized by a double-lobed phase trajectory, is shown in Fig. 4 (see below).

For other values of the drive frequency and amplitude, lying within this instability tongue, the motion of the pendulum can be more complicated. It turns out that this seemingly simple mechanical system exhibits an incredibly rich variety of nonlinear phenomena characterized by amazingly different regimes. Various regular and chaotic types of its dynamical behavior are described and partially explained in Ref. 11. Among them there are synchronized non-uniform unidirectional rotations in a full circle with a period that equals either the driving period or an integer multiple of this period, combined rotational and oscillatory motions, and motions synchronized (locked in phase) with oscillations of the pivot. Different competing modes can coexist at the same values of the driving amplitude and frequency. Which of these modes is eventually established, after the transient is over, depends on starting conditions; each limit cycle (attractor) has a definite basin of attraction. The motion of the pendulum can be irregular (chaotic), characterized by a strange attractor. Many interesting and even counterintuitive regimes are included as predefined examples in the simulation program^{10} “Pendulum with the vertically driven pivot.” Some of the limit cycles for the nonlinear pendulum with damping, which are easily reproduced in the simulations, help us in understanding the simple spectral composition of periodic oscillations that correspond to transition curves of the Ince–Strutt diagram for a frictionless system described by the linear Mathieu equation.

To calculate the critical driving amplitude (for a given drive frequency), which corresponds to the boundaries of instability “tongues,” we can restrict the theoretical analysis to infinitely small deviations of the pendulum from the equilibrium position. This allows us to replace $sin\u2009\phi $ by $\phi $ in the exact differential equation of the parametrically driven pendulum, Eq. (9), and omit the damping term. Thus Eq. (9) is reduced to the linear Mathieu equation, Eq. (2), in which the drive frequency *ω* is assumed as a unit of frequency

Periodic solutions of Eq. (11) that correspond to the desired boundaries of the instability interval, the interval that includes the principal parametric resonance, can be sought as a superposition of the fundamental harmonic whose frequency $\omega /2$ equals half the driving frequency, and a small admixture of the third harmonic component at the frequency $3\omega /2$,

This simple spectral composition of stationary oscillations at the transition boundaries (characterized by negligible contributions of higher harmonics) is clearly seen from the graphs in Figs. 3 and 4 (see below), obtained in numerical simulations.

Substituting the trial function $\phi (t)$, Eq. (12), into the Mathieu equation, Eq. (11), we find, from the requirement that the function (12) must satisfy Eq. (11), two separate (not coupled) systems of homogeneous equations for coefficients *C*_{1}, *C*_{3} and *S*_{1}, *S*_{3}, respectively.

The lower-frequency boundary of the principal parametric resonance (first instability tongue, transition curve 3 of the Ince–Strutt diagram, Fig. 1) corresponds to the cosine-like part of a non-trivial solution to Eq. (11), $\phi (t)=C1\u2009cos(t/2)+\u2009C3\u2009cos(3t/2)$. This periodic solution exists, if, for a given driving frequency *ω* (for a given value of parameter $k=(\omega 0/\omega )2=g/(l\omega 2)$ in this equation), the normalized amplitude of the pivot $m=a/l$ satisfies the following condition:

Periodic oscillations at this boundary are illustrated in Fig. 3. The graphs of $\phi (t)$ and $\phi \u0307(t)$ with their harmonic components are plotted for the drive frequency $\omega =1.529\omega 0$ and normalized amplitude $a/l=0.4$ (parameter *k* = 0.428 and *m* = 0.400). These values belong to the boundary 3 of the Ince–Strutt diagram defined in Eq. (13). In a linear system governed by the Mathieu equation, such cosine-like oscillations are excited at an arbitrary non-zero initial displacement and zero initial velocity. The relative contribution of the 3rd harmonic in this case is $C3/C1=\u22121+(4k\u22121)/(2m)=\u22120.112$.

Equation (13) gives an approximate analytical expression for the transition curve 3 of the Ince–Strutt diagram. A more precise expression for this curve 3 can be obtained with the trial function that includes also a (tiny) contribution of the fifth harmonic C_{5} cos(5*t*/2), see Sec. V.

Similarly, the higher-frequency boundary of the principal parametric resonance (curve 2 of the Ince–Strutt diagram, Fig. 1) corresponds to the sine-like part, $\phi (t)=S1\u2009sin(t/2)+\u2009S3\u2009sin(3t/2)$, of a non-trivial solution to Eq. (11). This solution exists, if, for a given driving frequency *ω* (for a given value of parameter $k\u22641/4$), the normalized amplitude of the pivot $m=a/l$ satisfies the following condition:

The relative contribution of the 3rd harmonic to periodic oscillations at this boundary is $S3/S1=1\u2212(1\u22124k)/(2m)$. Typical time-dependent graphs of the angular deflection $\phi (t)$ and of the angular velocity $\phi \u0307(t)$ (with the graphs of their harmonic components), and the double-lobed phase diagram for period-2 regular oscillations on the transition curve 2 are shown in Fig. 4. Actually, Fig. 4 corresponds to a stationary regime (limit cycle) of the nonlinear parametrically excited pendulum, described in Eq. (9) with a weak damping (*Q* = 100) in the absence of gravity (*g* = 0) and a slightly greater driving amplitude ($a/l=0.456$) than that for the boundary 2 of the Ince–Strutt diagram ($a/l=0.454$). In a linear system, such sine-like oscillations are excited at an arbitrary non-zero initial velocity and zero initial displacement.

Curves 2 and 3 that form the tongue of the principal parametric resonance in the Ince–Strutt diagram emanate at *m* = 0 from value *k* = 1/4 ($\omega =2\omega 0$). Rising from this *k*-value, curve 2 extends to smaller *k*-values (to higher excitation frequencies), intersects the ordinate axis *k* = 0 at $m=3(13\u22123)/4=0.454$, and continues further to negative *k*-values, that is, to negative values of *g*.

As we have already mentioned, we can treat the situation with negative values of parameter *g* in Eq. (9) as a field of “gravity” directed upward, which is exerted on the ordinary (hanging down) pendulum. In other words, we can treat the free-fall acceleration *g* in Eq. (9) as a control parameter whose variation is physically equivalent to the variation of the gravitational force. When this control parameter *g* is diminished through zero to negative values, the gravitational torque in Eq. (9) also reduces down to zero and then changes its sign to the opposite. Such a “gravity” tends to bring the pendulum into the inverted position $\phi =\pi $, destabilizing the lower equilibrium position $\phi =0$ of the unforced pendulum, and making upper position $\phi =\pi $ stable: at *g* < 0 the inverted position in Eq. (9) is equivalent to the lower position at positive *g*. That is, the region of negative *k*-values in the Ince–Strutt diagram corresponds formally to the stability of the inverted pendulum with vertically oscillating pivot in ordinary (directed downward) gravitational field.

The effect of the inverted pendulum dynamical stabilization by forced vertical vibration of the pivot (“Kapitza's pendulum”) is described by curve 1 (in the region *k* < 0) of the Ince–Strutt diagram, which shows the minimal amplitude *m* of the pivot that provides stability at given forcing frequency (see Sec. V below). If the amplitude exceeds the value given by curve 2, the (dynamically stabilized) inverted position becomes unstable through excitation of period-2 oscillations, which were described for the first time and experimentally demonstrated (the so-called “flutter” mode) in Refs. 12 and 13. Hence, curve 2 in the region of negative *k*-values corresponds to the upper boundary of the inverted pendulum dynamical stabilization. The “flutter” oscillations in the inverted pendulum and their relationships with the ordinary parametric resonance of the hanging down pendulum are discussed in detail in Ref. 14.

Intersection of curve 2 with the ordinate axis *k* = 0 occurs at $m=3(13\u22123)/4=0.454$. This value gives (for zero gravity) the amplitude of the pivot, which corresponds to transition from dynamically stabilized equilibrium of the pendulum to instability through excitation of period-2 oscillations like those at principal parametric resonance of the hanging pendulum, or like the “flutter” oscillations of the inverted pendulum. At *k* = 0 and *m* = 0.454 the ratio $S3/S1$ of amplitudes of the third harmonic to the fundamental one equals $\u2212(13\u22123)/6=\u22120.101$.

These values *m* = 0.454 and $S3/S1=\u22120.101$ agree well with the simulation experiment in conditions of the absence of gravity and very small angular excursion of the pendulum. When the normalized amplitude of the pivot $m=a/l$ exceeds the critical value $amin/l=0.454$, the angular excursion of the period-2 “flutter” oscillation (amplitude *S*_{1} of the fundamental harmonic) increases in proportion to the square root of this excess: $S1\u223ca\u2212amin\u2009$. This dependence follows from the nonlinear differential equation of the pendulum, Eq. (9), if the $sin\u2009\phi $ term in it is expanded as $\phi \u2212\phi 3/6$. This dependence also agrees well with the simulation experiment.

## IV. THE SECOND PARAMETRIC RESONANCE AND CURVES 4–5 OF THE INCE–STRUTT DIAGRAM

The second “tongue” of parametric instability in the absence of friction emanates at *m* = 0 from *k* = 1 ($\omega =\omega 0$). Stationary oscillations at the boundaries of this instability region are characterized by a dominant fundamental spectral component whose frequency *ω* equals the drive frequency, and several higher harmonics whose amplitudes rapidly diminish with their number. To find the relationship between the driving frequency and the amplitude of the pivot, which corresponds to existence of such periodic oscillations, we may substitute to Eq. (11) the trial function $\phi (t)$ in the following form:

The contributions of separate harmonic components (that is, amplitudes *C _{i}* and

*S*) can be found from the requirement that $\phi (t)$ in Eq. (15) satisfies Eq. (11). In the absence of friction, this requirement yields two independent systems of equations for

_{i}*C*and

_{i}*S*. This means that the periodic oscillations at one boundary of the instability tongue are of the cosine-type, while at the other of the sine-type. Equations for the corresponding bounding curves are given by the conditions of existence of non-trivial solutions to these homogeneous systems of equations.

_{i}For the sine-type boundary (curve 4 of the Ince–Strutt diagram) typical oscillations are shown in Fig. 5 (initial values $\phi (0)=0,\u2009\phi \u0307(0)\u22600$). If we omit the tiny contribution of the fourth harmonic component, $S4\u2009sin(4t)$, in the trial function φ(*t*) of Eq. (15), the equation is

If the contribution of the fourth harmonic component, $S4\u2009sin(4t)$, is included into the trial function $\phi (t)$ of Eq. (15), a more precise analytical expression for this curve can be obtained

The curves 4 plotted with the help of Eq. (16) or Eq. (17) are visually almost indistinguishable, so that for most applications of the Ince–Strutt diagram in physical problems distinctions between expressions (16) and (17) for curve 4 are practically insignificant.

For the cosine-type boundary (curve 5) typical oscillations are shown in Fig. 6 (initial values $\phi (0)\u22600$, $\phi \u0307(0)=0$). If we omit the tiny contribution of the third harmonic component, $C3\u2009cos(3t)$, in the trial function $\phi (t)$ of Eq. (15), the equation is

The same expression for this curve can be obtained otherwise as the limit $n\u2192\u221e$ in the general formula for the subharmonic resonance of order *n* that relates the drive amplitude $m(n,\u2009k)$ with the *k*-value (with the drive frequency), which gives the conditions of this resonance (see Sec. V for details).

If the contribution of the third harmonic component $C3\u2009cos(3t)$ is included in the trial function $\phi (t)$ of Eq. (15), a more precise analytical expression for this curve can be obtained

## V. SUBHARMONIC RESONANCES AND DYNAMICAL STABILITY OF THE INVERTED PENDULUM: CURVE 1

A physically transparent way to obtain the instability tongues for the pendulum whose pivot is forced to oscillate at a high frequency is based on the relationship between “slow” oscillations about the hanging down or the dynamically stabilized inverted position, and stationary periodic regimes of so-called subharmonic oscillations (see Ref. 15). This approach also allows us to find a more exact and enhanced criterion for dynamical stabilization of the inverted pendulum (see Ref. 16), which is valid over a wider region of the system parameters in comparison with a conventional criterion (see, for example, Refs. 9 and 14).

Within ranges of stability, the pendulum with damping, instead of gradually approaching the equilibrium position (either dynamically stabilized inverted position or ordinary downward position) by the process of damped slow oscillations, can be trapped into a *n*-periodic limit cycle, if initial conditions lie within a certain basin of attraction of this limit cycle. By virtue of the nonlinear effect of “phase locking” (synchronization of oscillations of the pendulum with those of the pivot, characterized by the certain relationship between their phases), the pendulum is regularly fed by additional energy to compensate for frictional losses.

The phase trajectory of such oscillations exactly repeats itself after *n* periods *T* of excitation are finished. One period of these non-damping oscillations of the dissipative pendulum equals exactly an integer number of cycles (*n*) of the pivot vibrations. The frequency $\omega /n$ of the principal harmonic equals $1/n$ of the excitation frequency *ω*. This allows us to call this phenomenon a subharmonic resonance of order *n*.

For the inverted pendulum with a vibrating pivot, periodic oscillations of this type were first described by Acheson and Mullin,^{17,18} who called them “multiple-nodding” oscillations. Computer simulations show that the pendulum motion in this regime has the appearance of some kind of original dance. Similar “dancing” oscillations can be executed also (at appropriate values of the driving parameters) about the ordinary equilibrium position (see Ref. 14). Figure 7 illustrates such “quadruple-nodding” regular oscillations about the lower equilibrium position. “Multiple-nodding” oscillations can occur also in the absence of gravity about any of the two equivalent dynamically stabilized equilibrium positions, see Ref. 15.

Understanding the spectrum of periodic oscillations occurring at subharmonic resonances will help us in finding relationships between the frequency and amplitude of the drive, at which such resonances take place. These relationships give the desired analytical expressions for some transitional curves of the Ince–Strutt diagram. In particular, the boundaries of the instability tongue of the principal parametric resonance (curves 2 and 3) correspond to conditions of the 2nd-order subharmonic resonance. Next we show that the lower-amplitude boundary of the dynamically induced stability of the inverted pendulum (curve 1) corresponds to the condition of *n*-order subharmonic resonance in the limit of infinitely large *n* ($n\u2192\u221e$). In other words, curve 1 corresponds to small-amplitude “slow” oscillations of the inverted pendulum about the upward vertical, whose period is indefinitely long.

To understand the spectrum of periodic subharmonic oscillations of an arbitrary order, next we consider a certain example of subharmonic resonances, namely, the resonance of the 6th order. The spectral composition of such oscillations is shown in Fig. 8. One period of the pendulum motion covers exactly six cycles of the pivot oscillations. Poincaré sections on the phase orbit mark the mechanical state ($\phi ,\u2009\phi \u0307$) of the pendulum once during each cycle of the excitation. The set of these sections (Poincaré map) consists of six fixed points in the phase plane visited by the representing point in a definite sequence. The fundamental harmonic component whose period equals six driving periods dominates the spectrum. We may treat this component of the spectrum as a subharmonic (as an “undertone”) of the driving oscillation.

Figure 7 and especially Fig. 8 show clearly that the spectrum of such *n*-periodic oscillations of a small amplitude consists of the principal harmonic $A1\u2009sin[(\omega /n)t]$ at the frequency $\omega /n$, two harmonics of order *n* − 1 and *n* + 1 (frequencies $(n\u22121)\omega /n$ and $(n+1)\omega /n$) with almost equal amplitudes, and a tiny admixture of two higher harmonics of order $2n\u22121$ and $2n+1$. Therefore, the time-dependence of $\phi (t)$ for the *n*-order subharmonic resonance can be sought in the following form:

Further on we assume the driving frequency *ω* as a unit of frequency (*ω* = 1), and substitute the trial function $\phi (t)$ of Eq. (20) into Eq. (9). Expanding the products of trigonometric functions, we obtain the following system of approximate equations for the contributions of the fundamental harmonic component and of the higher harmonics (for coefficients *A*_{1}, $An\u22121$ and $An+1,\u2009A2n\u22121$ and $A2n+1$):

The homogeneous system (21) has a nontrivial solution if its determinant equals zero. This condition yields an equation for the corresponding threshold (minimal) normalized driving amplitude $mmin=amin/l$ at which a certain *n*-periodic mode $\phi (t)$, given by expression (20), can exist. (We note that the same quadratic equation for $mmin2$ follows from the requirement of existence of non-zero solutions for the amplitudes of harmonic components, if in the trial function $\phi (t)$ of Eq. (20), all sine functions are replaced by cosine functions.) Then, after substituting this critical driving amplitude $mmin$ into the system (21), fractional amplitudes $An\u22121/A1,\u2009An+1/A1,\u2009A2n\u22121/A1$, and $A2n+1/A1$ of high harmonics for periodic oscillations of a given order *n* can be found as the solutions (for given values of *n* and *k*) to the homogeneous system of Eq. (21).

If we ignore in $\phi (t)$ the contribution of higher harmonics whose frequencies are $(2n\u22121)\omega /n$ and $(2n+1)\omega /n$, that is, assume $A2n\u22121$ and $A2n+1$ to be zero, system (21) simplifies considerably. The corresponding approximate solution can be found in Ref. 15. For the full system (21), the final expressions for $mmin$ and for the amplitudes of harmonics in Eq. (20) are too bulky to be cited here. We have used them in Fig. 9 for plotting the curves of $mmin$ as functions of $k=(\omega 0/\omega )2$ (inverse normalized driving frequency squared) corresponding to subharmonic oscillations of several orders *n* (namely, for *n* = 2, *n* = 4, and *n* = 6).

As a numerical example of such calculations, the following values of relative amplitudes $An\u22121=A5,\u2009An+1=A7,A2n\u22121=A11,\u2009A2n+1=A13$ of higher harmonics correspond to subharmonic oscillations of order *n* = 6 (shown in panel (a) of Fig. 10, see below) at *k* = 0, $m=0.2256\u2009l$:

These amplitudes refer to small oscillations of the pendulum in the absence of gravity, which are described by function $\phi (t)$ of Eq. (20), with all harmonics of the sine type. The phases of harmonics depend on initial conditions, and will be like those in Eq. (20) or Eq. (22), if the pendulum is excited by a small initial push from the equilibrium position ($\phi (0)=0$, $\phi \u0307(0)\u22600$). In Fig. 8, the relative contributions of harmonics satisfy Eq. (22), but their phases are different, because initially both $\phi (0)\u22600$ and $\phi \u0307(0)\u22600$ (mixed sine- and cosine-type oscillations).

The curves of subharmonic resonances of the order *n* = 2 give both boundaries of the instability tongue, which corresponds to the principal parametric resonance (curves 1–2 of the Ince–Strutt diagram). Analytical expressions (not cited here), which can be obtained from system (21) for these curves, are more complicated than those given in Eqs. (13) and (14). However, the difference in numerical values is very small. At least, the corresponding curves of the Ince–Strutt diagram plotted with the help of both analytical expressions visually coincide.

For subharmonic resonances of higher orders (*n* > 2), the sine- and cosine-like boundaries merge: the “tongues” of these resonances have zero width, and each of them is represented in the Ince–Strutt diagram by a single curve (see Fig. 9). This means that if the frequency and amplitude of the drive are chosen so that they correspond to subharmonic oscillations of a certain order *n* (lie on the relevant curve), such *n*-periodic stationary regime will be observed at arbitrary initial conditions. (The only requirement is that these initial conditions must provide oscillations of a small amplitude.) In particular, sine-like *n*-periodic stationary oscillations are excited at $\phi (0)=0,\u2009\phi \u0307(0)\u22600$, while cosine-like – at $\phi (0)\u22600,\u2009\phi \u0307(0)=0$.

For oscillations of a finite amplitude, when nonlinear properties of the pendulum become essential and its motion obeys the exact differential equation, Eq. (9), the above-described subharmonic oscillations transform into *n*-periodic limit cycles with a similar spectrum. The pendulum with damping eventually becomes trapped in this regular periodic motion if the initial conditions lie within its basin of attraction.

We note that the curves in the Ince–Strutt diagram, which correspond to periodic oscillations of a certain order *n*, whose spectrum consists of the principal harmonic at the frequency $\omega /n$ and higher harmonics of orders *n* − 1 and *n* + 1 (frequencies $(n\u22121)\omega /n,\u2009(n+1)\omega /n$, and a tiny admixture of harmonics of order $2n\u22121$ and $2n+1$), have branches embedded in stable regions between curves 1–2, between curves 3–4, and also in other stable regions. In particular, regular oscillations of the period 6*T* (*n* = 6) occur, if the frequency and amplitude of the pivot (*k* and *m* values) lie on the curves labeled *a*, *b,* and *c* in Fig. 9. Graphs of the time dependence and the spectra of these oscillations are shown in panels (a), (b) and (c) of Fig. 10, respectively.

For curve *a*, the fundamental harmonic whose frequency is $\omega /6$ (an “undertone” of the drive frequency *ω*) dominates the spectrum. In our physical model, this harmonic describes the smooth (slow) motion of the pendulum in the effective potential well, which is created by fast oscillations of the pivot.^{14} The smooth motion is distorted by a small superposition of high harmonics at frequencies $5\omega /6$ and $7\omega /6$.

In the spectrum of period-6 oscillations that correspond to curve *b* of Fig. 9, the contribution of harmonic component at the frequency $(5/6)\omega $ is the greatest. This means that during one period of such oscillations, which equals exactly 6 drive periods *T*, the pendulum executes 5 cycles of its motion.

Similarly, for curve *b* of Fig. 9, the harmonic component at frequency $(7/6)\omega $ has the greatest amplitude. Therefore, the period 6*T* of these oscillations spans over 7 cycles of the pendulum's motion.

The analytical expression for curve 1 of the Ince–Strutt diagram can be obtained from the general formula for the subharmonic resonances (too complicated to be cited here), which relates the frequency *ω* [or parameter $k=g/(l\omega 2)$] with the normalized amplitude $m=a/l$ of the pivot. This can be done (using computer algebra) by the transition to the limiting case $n\u2192\u221e$ in the formula mentioned. Indeed, the lower boundary of dynamical stabilization of the inverted pendulum corresponds to the transition from unstable equilibrium to the dynamically stabilized one through neutral equilibrium, in which the upside-down pendulum moves indefinitely slowly in the vicinity of the inverted position, that is, swings about the upward vertical with an infinitely long period. This regime of infinitesimal oscillations at almost zero frequency ($\omega \u21920$) can be treated as a subharmonic oscillation of order *n* at the frequency $\omega /n$ in the limit $n\u2192\u221e$. Transition to this limit in the general formula yields a simple analytical expression, Eq. (3), for curve 1 of the Ince–Strutt diagram, Fig. 1.

## VI. TRANSITION CURVES OF THE INCE–STRUTT DIAGRAM IN THE CASE OF DAMPED MATHIEU EQUATION

Our approach to finding analytical expressions for the transition curves in the Ince-Strutt diagram by considering the corresponding periodic regimes of a parametrically driven pendulum can also be used for the physical system described by the Mathieu equation with a damping term, namely, as in Eq. (9). The boundaries of the periodic modes correspond to oscillations with infinitesimal amplitudes ($\phi \u226a1$), so that when finding these boundaries we can use the linearized equation with damping

To characterize the intensity of friction, we have introduced here (instead of *γ*) another (dimensionless) parameter *β*,

Here, $Q=\omega 0/2\gamma $ is the dimensionless quality factor, commonly used to characterize the intensity of viscous friction.

To find the instability tongue of the principal parametric resonance, which occurs at excitation frequencies *ω* approximately twice the natural frequency *ω*_{0} (resonance of order *n* = 1), we substitute the trial function $\phi (t)$, Eq. (12), into the Mathieu equation with damping, Eq. (23). From the requirement that the function (12) must satisfy Eq. (23), we find the system of homogeneous equations for coefficients *C*_{1}, *C*_{3} and *S*_{1}, *S*_{3},

A non-trivial solution to this system exists if its determinant equals zero. This condition can be obtained with the help of computer algebra and yields an equation (not cited here) that defines the dependence of driving amplitude $m=a/l$ on *k* (on frequency *ω* of the excitation) for periodic motions that correspond to the desired boundaries of the instability tongue. Contrary to the frictionless case, the tongue in the Ince–Strutt diagram (Fig. 11) does not touch the abscissa axis, because in the system with friction parametric resonance is possible if the normalized amplitude of excitation $m=a/l$ exceeds some threshold (minimal) value $mmin$, which is inversely proportional to the quality factor *Q*. For the principal resonance, the threshold amplitude is approximately evaluated as $mmin\u22481/(2Q)$. For example, at exact tuning to the principal resonance (*k* = 1/4, that is, $\omega =2\omega 0$) the normalized axis amplitude $m=a/l$ of the pivot at *Q* = 5.00 must be greater than 0.100. The above-mentioned analytical equation for the tongue boundary yields the threshold *m*-value at the exact resonance with a somewhat greater precision, namely, 0.1002. However, according to this same analytical equation, the lowest point of this instability tongue is slightly shifted from the exact resonance: it is located at *k* = 0.244 ($\omega =2.025\omega 0$) and corresponds to the absolute threshold value $mmin=0.099$ of the drive amplitude for the system with *Q* = 5.00.

To find the second tongue of parametric instability ($\omega \u2248\omega 0$) in the presence of friction, we may substitute into Eq. (23) the trial function $\phi (t)$ given in Eq. (15). This yields the following system of homogeneous equations for contributions *C _{i}* and

*S*of certain harmonics into the stationary period-1 oscillations, corresponding to the boundary of the second tongue:

_{i}From the condition of existence of a non-trivial solution to this system, we get an equation (not cited here) for the desired boundary of the second instability tongue. Figure 11 shows the corresponding transition curve, plotted according to this equation for the system with *Q* = 5.0. The bottom point of this tongue is located at *k* = 1.106 ($\omega =0.951\omega 0$) and *m* = 0.999. That is, the threshold amplitude *a* for excitation of the second parametric resonance in a system with friction, characterized by the quality factor *Q* = 5.0, is about ten times greater than for the principal resonance: it equals the pendulum length *l*. The graphs of periodic oscillations at the threshold of second-order parametric excitation (regime of parametric regeneration) are shown in Fig. 12. For the second resonance, the threshold (minimal) value $mmin$ is approximately evaluated to be inversely proportional to the square root of the quality factor *Q*: $mmin\u22482.32/Q$.

## VII. CONCLUSIONS

Many problems in physics and engineering that involve the Mathieu equation reduce to the determination of conditions, in which solutions to this equation either remain bounded or grow indefinitely in the course of time. The answer is given by the Ince–Strutt diagram, which shows the transition curves in the parameter plane between different regimes. Stability charts of the Mathieu equation published in the literature are usually calculated using various numerical methods.

However, it is possible to obtain for the instability boundaries fairly simple finite-form analytical expressions, which can be helpful in the investigation of various physical systems described by the Mathieu equation. In this paper, we have shown how the derivation of such expressions can be done by referring to the physical sense of transition curves in the Ince–Strutt diagram for a familiar and well-studied mechanical system—the rigid planar pendulum whose pivot is forced to oscillate along the vertical line.

This parametrically excited physical system is described by a nonlinear equation with damping, which is a generalization of the Mathieu equation. For such a nonlinear damped pendulum, various periodic regimes can be easily reproduced and observed in computer simulations. Stationary oscillations of the pendulum in such periodic regimes are characterized by rather simple spectral compositions, which are easily determined in simulations. Knowledge of the spectrum allows us to find (using a modification of the harmonic balance method) the desired analytical expressions for the transition curves in a system described by a simpler differential equation (linear Mathieu equation without damping), for which it is not so easy to reproduce periodic regimes in simulations. The desired transition curves in a linear system correspond to periodic oscillations, which include the same harmonic components as similar periodic regimes (limit cycles), easily observed for the nonlinear pendulum with damping. It thus turns out that simulations of a somewhat more complex nonlinear system can help us in investigations of linear systems described by the Mathieu equation.

## References

*Nonlinear Oscillations*, Computer simulations (