## Abstract

This article investigates the equilibria and stability of a pendulum when the support has a prescribed motion defined by an elliptic function. Stability charts are generated in the parameter plane for different values of the elliptic function modulus. Numerical integration and Floquet theory are used to generate stability charts that are later obtained through harmonic balance analysis. It is shown that the size and location of the instability tongues is directly linked to the elliptic function modulus. Comparisons are also made between the stability charts of Mathieu's equation and those of the pendulum when the prescribed motion is defined by an elliptic function.

## 1. Introduction

A physical system undergoes a parametric forcing if one of its parameters is modulated periodically with time. A familiar example of parametric excitation is given by a pendulum whose support is given a prescribed motion in the vertical direction. It is known that when the frequency of a sinusoidal driver is twice the natural frequency of the pendulum, the downward equilibrium position becomes unstable, and the oscillation amplitude increases exponentially. This phenomenon is called parametric resonance. To investigate this phenomenon and thus the stability of the pendulum's downward equilibrium, one typically linearizes the equation of motion about the equilibrium point, which gives the following differential equation: 1.1

The same equation holds for the stability of the upward equilibrium position, although *δ* becomes negative. Equation (1.1) is known as Mathieu's equation and it was first discussed by Mathieu (1868) when he researched the vibration of an elliptic membrane. The stability or boundedness of solutions of Mathieu's equation (1.1) as a function of parameter values (*ϵ*,*δ*) has been extensively examined and characterized in terms of transition curves. These curves consist of one-dimensional curves in the parameter plane (*ϵ*,*δ*) and they separate regions of stability from regions of instability. In order to study the stability of equation (1.1), Hill (1886) presented a method called Hill's infinite determinant method. Other techniques have been used to investigate the stability of equation (1.1), such as Lyapunov–Floquet transformation, and perturbative Hamiltonian normal forms (Butcher & Sinha 1995), integral of energy and numerical integration (Jazar 2004), and other perturbation methods that involve a small parameter (Stoker 1950; Nayfeh & Mook 1979).

Many authors have investigated the stability chart of systems that involve Mathieu's equation. Rand *et al.* (2010) investigated the effect of fractional damping on the transition curves in Mathieu's equation (1.1). It was shown that the shape and location of the *n*=1 transition curve can be changed by changing the value of the order of the fractional derivative. In earlier work, a linear damping term was considered in Taylor & Narendra (1969), Gunderson *et al.* (1974) and Richards (1976), and it was shown that adding a damping term increases the area of stability in parameter space. Insperger & Stépán (2002) studied equation (1.1) with a delay term and constructed a closed-form three-dimensional stability chart for the delayed Mathieu equation. They found analytically that the boundary curves in the plane (*δ*,*b*), where *b* is the amplitude of the delay term, are lines for any *ϵ*. These boundary curves delimit triangles of stability that diminish as *ϵ* increases. Rand and co-workers (Zounes & Rand 1998; Rand *et al.* 2003; Rand & Morrison 2005) considered equation (1.1) with a second parametric forcing function and investigated the transition curves in the quasi-periodic Mathieu equation. It was shown that each of the instability regions corresponds to a distinct order of resonance between the slow flow motion of equation (1.1) and the frequency of drift between the two drivers. It was also shown that the instability region near 2 : 1 : 1 resonance is much thicker than the comparable region near 2 : 2 : 1 resonance.

In equation (1.1), the forcing function is of trig type; however, systems may be driven by forcing functions of elliptic type. For instance, Sanjuán (1998) studied a nonlinear pendulum parametrically forced with a Jacobi elliptic function, and showed that the periodicity of a solution can be switched by changing the wave form and the periodicity of the driving force. Cveticanin (2000) considered a strongly nonlinear system parametrically excited with a Jacobi elliptic function, where the elliptic-Krylov–Bogolubov method was used to approximate regions of unbounded solution. Elías-Zúñiga (2004) investigated the effect of applying an axial force of elliptic type to a hinged beam–column system by studying the exact solution of Lamé's equation. For a detailed treatment of Lamé's equation, see Arscott (1964). In a similar problem to Elías-Zúñiga, Artim & Aydin (2010) considered an axial loading of elliptic type, but with a force that combines three Jacobi elliptic functions.

This article investigates the qualitative and quantitative changes in the stability behaviour when the trigonometric term in equation (1.1) is replaced by the elliptic function cn(*t*,*m*). This work is motivated by the curiosity of the authors and the important role that elliptically forced systems could play in the field of nonlinear systems. For example, previous work has shown that elliptic forcing could be used as a way to suppress chaotic behaviour. This suppression is achieved by changing the shape of the periodic forcing, that is changing the modulus of the elliptic function (Chacońn & Dõńaz Bejarano 1993; Rodelsperger *et al.* 1995). Although much work has been done on systems that are parametrically excited with a Jacobi elliptic function, to the authors' knowledge, this is the first work where an approximation of the transition curves for a linear system excited parametrically with a Jacobi function is given. The approximations of the transition curves presented in this work could help to understand the dynamics of nonlinear systems with elliptic parametric forcing, since these curves provide global information on the stability of all equilibrium points in the case of a nonlinear system.

The remainder of this paper is organized as follows. Section 2 describes the model for a planar pendulum when the pivot point has a prescribed motion in the form of an elliptic function. Stability boundaries are then obtained from numerical studies of equation (2.3). Analytical studies are then performed to obtain approximate expressions for the transition curves. Trends in the stability behaviour are then discussed in §5 prior to offering concluding remarks in §6.

## 2. Parametrically excited pendulum

Consider the planar pendulum shown in figure 1 whose support has a prescribed motion defined by an elliptic function. The governing differential equation is
2.1where cn(*t*,*m*) is the cosine Jacobi elliptic function (Abramowitz & Stegun 1972) and *m* is its squared modulus, 0≤*m*<1. The amplitude and frequency of excitation and time *t*, have been rescaled,
2.2In order to investigate the stability of one of the equilibrium solutions *x*=0 or *x*=*π* of the pendulum, we linearize equation (2.1) about the desired equilibrium, giving an equation of the form
2.3The parameter *δ* is positive for the *x*=0 equilibrium position, and negative for *x*=*π*. The stability of equation (2.3) is symmetric with respect to *ϵ*=0. Our goal in this article is to investigate the transition curves in equation (2.3). These curves separate regions of stability from regions of instability. Here, stability is defined by all solutions being bounded, and instability by the existence of an unbounded solution. We note that for *m*=0 , equation (2.3) reduces to Mathieu's equation (1.1). Therefore, figure 2 displays the transition curves in equation (2.3) when *m*=0. Equation (2.3) will be referred to as the elliptic Mathieu equation. In §3, we present the stability charts in equation (2.3) obtained numerically.

## 3. Numerical integration

In this section, we present stability charts of the elliptic Mathieu equation by numerically integrating equation (2.3). In order to numerically investigate the transition curves in equation (2.3), we begin by replacing *t* with a new time variable *τ*=am(*t*,*m*) (Magnus & Winkler 1966). This transformation of time, which was also used in Recktenwald & Rand (2007), turns out to convert the Jacobi elliptic function in equation (2.3) to trig functions resulting in the following system:
3.1where
3.2thus
3.3Since equation (3.1) has periodic coefficients, this allows us to use numerical integration in conjunction with Floquet theory (Rand 2005). The stability of the system is found by creating a fundamental solution matrix. The differential equation (3.1) is numerically integrated for exactly one period of the forcing function, *T*=2 *π*, using two sets of initial conditions,
3.4The results are combined into a fundamental solution matrix *C* whose eigenvalues determine stability,
3.5If, for a given point in parameter space, either eigenvalue of *C* has absolute value greater than unity, then that point is said to be unstable (and stable otherwise). Figure 3 shows the stability charts of equation (3.1), consequently equation (2.3), for two different values of the modulus *m*=(0.5,0.95). From figure 3, we note that the number of the instability tongues (white regions) increases as the modulus *m* increases, and the size of these instability tongues changes as *m* is increased. Also, the measure of the stable region that lies in the *δ*<0 half-plane decreases as *m* increases. In order to analytically analyse the stability charts obtained numerically, we look for approximations of the transition curves by using the method of harmonic balance in §4.

## 4. Harmonic balance

In order to obtain an analytical approximation for the transition curves in equation (2.3), we apply the method of harmonic balance to equation (3.1). The reader can refer to Rand (2005) for a detailed description of the method of harmonic balance. From Floquet theory, it is known that there exist periodic solutions of period *T* and 2*T* on the transition curves, where *T* is the period of the forcing function, *T*=2*π* in equation (3.1). Therefore, we look for solutions along the transition curves in the form of Fourier series
4.1This series represents a general periodic function with period 4*π* and, as a special case, functions with period 2*π*. Substituting equation (4.1) into equation (3.1), using Macsyma to perform a trigonometric reduction and collect terms gives four sets of algebraic equations on the coefficients *a*_{n} and *b*_{n} (even and odd *n*'s).

Figure 4 shows the transition curves for *n*=80, and for modulus values *m*=(0.5,0.95). The analytical results obtained using the method of harmonic balance (figure 4) agree with the numerical results (figure 3). It should be noted that if some of the tongues in figure 3 seem to be detached from the *δ*-axis, it is just a numerical artefact arising from the integration algorithm used by Matlab.

The algebraic equations relating *m*, *ϵ* and *δ* are too long to list here. However, in order to track the *n*=1 transition curves as *m* changes, we posit a truncated Fourier series,
4.2After substituting and collecting terms, we obtain the following approximate expression for the *n*=1 transition curve:
4.3If *m* in equation (4.3) is replaced by *m*=0, we obtain an approximate expression of the *n*=1 transition curve in Mathieu's equation (1.1). Figure 5*a* shows how the location of the *n*=1 transition curve on the *δ*-axis changes as *m* varies. As *m* is increased, the instability tongue is translated to the left-hand side from the value 1/4 on the *δ*-axis. We note that in figure 5*a*, we can only see the change of the location; however, the shape of the instability tongue also changes with *m* (figures 3 and 4). In order to track the change of the shape, more harmonics need to be included in the ansatz (4.2).

An approximate expression of the *n*=0 transition curve is obtained by looking for a solution in the following form:
4.4Proceeding as before, we obtain the following approximate expression of the *n*=0 transition curve:
4.5Figure 5*b* shows the *n*=0 transition curve. We note that the shape of the transition curve changes as *m* varies. As *m* is increased, the *n*=0 transition curve moves to the left-hand side from the transition curve corresponding to *m*=0.

An interesting region in the (*ϵ*,*δ*) plane is the stability region bounded by the *n*=0 and *n*=1 transition curves. From figures 3 and 4, we note that as *m* increases, the measure of the stable region that lies in the *δ*<0 half-plane decreases. This can be shown analytically by solving the *n*=0 transition curve, equation (4.5), for *δ* and subtracting it from the expression for the left *n*=1 transition curve, equation (4.3), which gives, after developing in the Taylor series the following approximation:
4.6Therefore, for a fixed *ϵ*, the width of the stable region diminishes as *m* tends to 1.

Also, from figures 3*a* and 4*a*, we note that for *m*=0.5, the *n*th resonance tongue has a tip sharper than the (*n*−1)th one, where *n*=(1,2,3). However, for *m*=0.95, figures 3*b* and 4*b* show that the tongue tips become narrower as the resonance number *n* increases, except for the *n*=3 resonance tongue whose tip is broader than the one at *n*=2.

For example, at *n*=1, the boundary curves are intersecting at positive angles at *ϵ*=0, equation (4.3), while at *n*=2, where the transition curves have the following approximate expressions:
4.7and
4.8the boundary curves are tangents of order 1, i.e. with a contact of order 2. Now at *n*=3, the expression of the transition curves is
4.9where *δ*_{3,i} are listed in appendix A. From this expression, equation (4.9), we note that the order of contact is 1 since the coefficient at order *ϵ* is different for the two boundary curves; when *m*=0, *δ*_{3,1}=0 and thus we recover the order of contact at *n*=3 for the case of the Mathieu equation, which is 3. Figure 6 shows the instability width to a leading order for the *n*=2 (dashed) and *n*=3 (solid) resonance tongues that are, respectively, given by
4.10and
4.11

From this figure, we remark that as *m* is increased, the *n*=3 resonance tongue tip gets broader than the *n*=2 resonance tongue tip (figures 3 and 4). It should be noted that even though the *n*=1 and *n*=3 resonance tongues have the same order of contact, the coefficient at the order *ϵ* for the *n*=3, equation (4.9), is smaller than the one for the *n*=1, equation (4.3); this is why the *n*=3 resonance tongue tip is narrower than the *n*=1 one.

Equations (4.7)–(4.9) were obtained by expressing the transition curves in the form of a power series expansion
4.12and then substituting it into the four determinants obtained by the harmonic balance method, and finally equating the coefficients of the same order of *ϵ* to obtain the *δ*_{i} (Rand 2005).

## 5. Discussions

From figures 3 and 4 (*m*=0.5,0.95), we note that the number of ‘resonance tongues’ increases as the modulus *m* increases in the range 0≤*m*<1. Also the size of these instability tongues changes as *m* increases. To understand how the location of the resonance tongues on the *δ*-axis changes as *m* varies, we use Floquet theory to equation (2.3) as it has a periodic forcing for 0≤*m*<1 with period *T*=4 *K*(*m*), where *K*(*m*) is the complete elliptic integral of the first kind (Abramowitz & Stegun 1972),
5.1Floquet theory tells us that equation (2.3) for 0≤*m*<1 has solutions of period 4 *K*(*m*) and 8 *K*(*m*) on its transition curves. Now when *ϵ*=0, equation (2.3) takes the form d^{2}*x*/d*t*^{2}+*δ* *x*=0 and has solutions of period . These will correspond to solutions of period 4 *K*(*m*) and 8 *K*(*m*) when , which gives
5.2For *m*=0, *K*(0)=*π*/2, equation (5.2) becomes *δ*=*n*^{2}/4, which are the already known resonant values of *δ* in Mathieu's equation (1.1) (Rand 2005). Now for *m*=0.95 and *n*=5, equation (5.2) gives *δ*=1.8232, see the rightmost resonance tongue in figure 4*b*. Thereby, varying *m* changes the resonant values of *δ* in equation (2.3), in particular, as *m* increases the number of the resonant *δ*'s increases in a specified range of *δ*. Figure 7 displays the change of the resonant *δ*'s as a function of *m*, equation (5.2). From this figure, we note that the resonant values of *δ* approach zero as *m* tends to 1; this explains the increasing number of resonance tongues in the *δ*–*ϵ* parameter plane.

On the other hand, the transition curve emanating from *n*=0 does not change location on the *δ*-axis, whereas it changes shape as the modulus *m* varies (figure 5*b*).

From the results obtained by studying the elliptic Mathieu equation (2.3), we note that the regions where the equilibrium solutions *x*=0 are unstable change size as the modulus *m* varies (figures 3 and 4). The parametric resonance condition also changes with the modulus *m*, for example, in order to actuate the first parametric resonance, *n*=1, the frequency of the driver should be
5.3where in equation (5.2) and where is the natural frequency of the pendulum near the *x*=0 equilibrium position. Therefore, as *m* is increased, the frequency of the driver *w* needs to be increased in order to destabilize the downward equilibrium position. We also note that in order to hit the instability tongue at the *n*=3 parametric resonance, the amplitude of the driver needs to be increased as the modulus *m* decreases, since the *n*=3 resonance tongue tip gets sharper as *m* decreases.

In the case of the *x*=*π* equilibrium solutions, the region where these solutions are stable, i.e. on the left of *δ*=0 in figures 3 and 4, decreases in size as *m* increases. Therefore, increasing *m* decreases the possibility of stabilizing the upward equilibrium position of the pendulum, equation (4.6).

## 6. Conclusion

In this paper, we have investigated the stability of the downward and upward equilibria of a pendulum whose support is vertically excited with a force of elliptic type. The elliptic Mathieu equation (2.3) was studied to deduce the stability of the equilibrium positions of the pendulum. The stability charts of the elliptic Mathieu equation were produced by using numerical integration in conjunction with Floquet theory. The transition curves separating regions of stability from regions of instability in equation (2.3) were obtained by using the method of harmonic balance. We found that by increasing the value of the modulus *m* of the elliptic forcing function, the number of the instability tongues that corresponds to the regions where the *x*=0 equilibrium is unstable increases, and the size of these tongues changes. We gave approximate expressions of the *n*=0,1,2,3 transition curves. We also gave the expression of the resonant values of *δ* from which the transition curves emanate. It was shown that the *n*=3 resonance tongue changes order of contact from 3 to 1 when *m*≠0, also as the modulus *m* increases, the *n*=3 resonance tongue tip becomes broader than the *n*=2 resonance tongue tip.

A particular interest was the comparison of equation (2.3) with equation (1.1). It was found that the size of the instability tongues in the elliptic Mathieu equation (2.3) is narrower than the ones in the Mathieu equation (1.1).

It was also shown that as *m* is increased, the frequency of the driver *w* needs to be increased in order to destabilize the downward equilibrium position. On the other hand, we found that in order to hit the instability tongue corresponding to the 3 : 1 subharmonic resonance, the amplitude of the driver needs to be increased as the modulus *m* decreases. Also, increasing *m* decreases the possibility of stabilizing the upward equilibrium position of the pendulum.

## Acknowledgements

S.S. thanks Richard Rand for his helpful comments and his support. B.M. gratefully acknowledges the support from US National Science Foundation CAREER Award (CMS-0348288).

## Appendix A

A1

- Received May 28, 2012.
- Accepted August 24, 2012.

- This journal is © 2012 The Royal Society