## Abstract

An inverted pendulum with asymmetric elastic restraints (e.g. a one-sided spring), when subjected to harmonic vertical base excitation, on linearizing trigonometric terms, is governed by an asymmetric Mathieu equation. This system is parametrically forced and strongly nonlinear (linearization for small motions is not possible). However, solutions are scaleable: if *x*(*t*) is a solution, then so is *αx*(*t*) for any real *α*>0. We numerically study the stability regions in the parameter plane of this system for a fixed degree of asymmetry in the elastic restraints. A Lyapunov-like exponent is defined and numerically evaluated to find these regions of stable and unstable behaviour. These numerics indicate that there are infinitely many possibilities of instabilities in this system that are missing in the usual or symmetric Mathieu equation. We find numerically that there are periodic solutions at the boundaries of stable regions in the parameter plane, analogous to the symmetric Mathieu equation. We compute and plot several of these solution branches, which provide a relatively simpler means of computing the stability transition curves of this system. We prove theoretically that such periodic solutions must exist on all stability boundaries. Our theoretical results apply to the asymmetric Hill's equation, of which the pendulum system is a special case. We demonstrate this with numerical studies of a more general asymmetric Mathieu equation.

## 1. Introduction

Consider an inverted pendulum with asymmetric elastic restraints in which the pivot is given a vertical periodic oscillation, as shown schematically in figure 1. Spring stiffnesses on the left and right are taken as and , respectively, with . The springs may be taken as straight, as shown, with a long free length so that they remain effectively horizontal even when their endpoints move up and down; they may alternatively, for simplicity, be thought of as torsion springs at the pivot point. Gravity may be incorporated in the spring stiffnesses and is not explicitly accounted for. There is no damping.

Let the angular displacement of the pendulum be *θ*(*t*). Taking the length of the pendulum to be *l* and the origin of the coordinate system at some point on the *θ*=0 line, the coordinates of the centre of mass of the pendulum arewhere *u*(*t*) is the imposed displacement of the pivot. Differentiating,The system kinetic energy iswhere *m* is the mass. We take the potential energy to be of the form (recall that gravity is not explicitly included)which is correct for torsion springs, and approximately correct for straight springs with long free lengths as mentioned earlier.

The equation of motion can now be routinely obtained using Lagrange's method. On linearizing trigonometric terms by substituting cos *θ*=1 and sin *θ*=*θ*, and letting *l*=1, *m*=1, we obtain

Assuming *u*=*ϵ* cos(*t*), writing *x* in place of *θ*, and noting thatwe obtain(1.1)We call equation (1.1) an asymmetric Mathieu equation; setting gives the usual Mathieu equation (Stoker 1950; Winkler 1966; Rand 2003). We mention in passing that the dynamics of a gear-pair system involving backlash can be modelled by a similar equation (Theodossiades & Natsiavas 2000). In this paper, we investigate some of the similarities and differences between the stability diagrams (in the *δ*–*ϵ* parameter plane) for the usual and asymmetric Mathieu equations. Our theoretical treatment allows parametric forcing of the term as well, which introduces two more free parameters as discussed later.

As is well known, the usual Mathieu equation has alternating stable and unstable regions in the parameter plane. On the transition curves (boundaries of the stability regions), there are periodic solutions with period either 2*π* or 4*π* (2:1 subharmonic motion) (Stoker 1950; Winkler 1966; Rand 2003). How are these related to the stability regions, transition curves, and periodic solutions of the asymmetric Mathieu equation?

Equation (1.1) is a homogeneous second-order differential equation with real periodic coefficients. It is also essentially nonlinear: linearization near *x*=0 is not possible. If *x*_{1}(*t*) and *x*_{2}(*t*) are two solutions, then *αx*_{1}(*t*)+*βx*_{2}(*t*) is not a solution in general. However, solutions are *scaleable*, i.e. if *x*(*t*) is a solution, then *αx*(*t*) is also a solution for any real *α*>0. This scaleability will be important in the subsequent analysis.

## 2. Unforced system

The unforced system (*ϵ*=0) is(2.1)The asymmetric potential energy of this unforced system and the corresponding asymmetric spring force are shown in figure 2, left and right, respectively.

All solutions of the unforced system are periodic, and have the same period. Each such periodic solution consists of two half-sinusoids, one for *x*>0 and one for *x*<0, with different amplitudes and time durations. The time period *T* and corresponding fundamental frequency *ω* are(2.2)

## 3. A Lyapunov-like exponent for scaleable systems

Since the fully nonlinear forced system seems analytically intractable, we resort to numerics. A few numerically obtained solutions for equation (1.1) are shown in figure 3. It is seen that solutions grow without bound (unstably) for some parameter values, but remain bounded for others. The unstable solutions grow exponentially, and we begin by characterizing their exponential growth rate.

A Lyapunov-like exponent is defined as(3.1)where the state vector ; denotes the Euclidean norm; the subscript *k* in denotes the time in number of forcing periods; and a superscript, either 0 or 2*π*, denotes the start or end of a forcing period, respectively.

The above Lyapunov-like exponent is calculated as follows. We numerically integrate equation (1.1) over an interval of 2*π* with random initial conditions that satisfy . At the end of this integration, we have the state . Initial conditions for the next cycle are then obtained by rescaling to unit norm:Integration over another period of forcing then gives . The above steps are repeated *N* times, for some large *N*. We discard the first several states in the calculation to get rid of transients, and take the first retained state as . *N* is taken large enough to obtain convergence, which in this case is relatively rapid presumably because solutions are not chaotic.

An exponent *σ*=0 implies the corresponding *δ*–*ϵ* point is stable, while *σ*>0 implies it is unstable. In numerics, *σ*=0 is not exactly observed, and a non-zero numerical tolerance or threshold is set by the analyst.

## 4. Numerical stability diagram

We first present the results of a numerical stability analysis of the asymmetric Mathieu equation on the *δ*–*ϵ* parameter plane for the arbitrarily chosen value of (we will show some more results for later in this same section). In these calculations, numerical integration of ODEs was done using Matlab's built-in ODE solver ‘ODE45’ with ‘event detection’ so as to ensure that every zero-crossing of *x* coincides with a mesh point (or sampling instant) in the time-discretized solution.

Figure 4 shows the Lyapunov-like exponent on a 500×500 grid covering *δ*∈[−0.1,1.3] and *ϵ*∈[0,1], for *N*=600. The calculation took several days, distributed among a few PCs. Figure 5 shows a closer and refined view of the boxed, bottom left, region of figure 4. Now *δ*∈[0,0.45] and *ϵ*∈[0,0.4], the grid is 450×400, and *N*=2400.

The figures show several instability regions. Limitations in numerics and graphics limit the number of visible regions; however, it seems that there are infinitely many such regions, successively narrower and more weakly unstable, nested between a few large and strongly unstable regions. These large and strongly unstable regions have a one-to-one correspondence with regions of instability for the usual Mathieu equation, as may be seen from the following. Narrower and finer resonance regions are associated with smaller numerical values of Lyapunov-like exponents: decreasing *N* or increasing the numerical value of the cutoff tolerance in the stability diagram will prevent us from detecting these regions.

For the usual or symmetric Mathieu equation (), there are regions of instability emanating from the *δ*-axis at *δ*_{sym}=0, 0.25, 1, ⋯. These correspond to simple rational relationships between the unforced time period, , and the forcing period, 2*π*.

For the asymmetric Mathieu equation with , the time period of the unforced system is, from equation (2.2),For the same rational relationships to hold so that the same resonances may occur, we require this time period to bear the same ratio to 2*π* as in the symmetric case.

This givesorThus, we expect corresponding instability regions for the case to emanate from the *δ*-axis at *δ*=0, 0.4201, 1.6805, ⋯. This is supported by the numerical results.

The asymmetric and symmetric Mathieu equations are identical on the *ϵ*-axis (with *δ*=0),

Therefore, stability intervals on the *ϵ*-axis for both equations (for any ) are the same. So, all the new instability regions of the asymmetric equation must have zero widths on the *ϵ*-axis. This is observed in the numerics. Similarly, the widths of the corresponding strong instability regions that occur for both the asymmetric as well as symmetric equations coincide on the *ϵ*-axis; this is outside the region covered by the numerical results, and cannot be seen in the figure.

We now consider other possible resonances. Let *r* be a simple rational number, such as (say) 2.5, 3 or 4. Consider a situation where the ratio between the unforced time period and the forcing period is *r*, i.e.(4.1)orThe *δ* values corresponding to *r*=2.5, 3 and 4 are then 0.2689, 0.1867 and 0.1050, respectively. There are, in fact, instability regions corresponding to these values on the *δ*-axis, as may be seen with a little faith from figure 5, or more convincingly from figure 7 in §5.

The qualitative results obtained above are not special for . Figure 6 shows the Lyapunov-like exponent for on a 500×500 grid covering *δ*∈[−0.1,1.5] and *ϵ*∈[0,1], for *N*=600. There is agreement with the results for .

From the numerical results of this section, it seems likely that all the narrow instability regions seen in figure 5 do in fact continue all the way to the *δ*-axis. This is verified numerically in §5, and supported further with theory later in the paper.

## 5. Periodic solutions on stability boundaries

As is well known, for the usual or symmetric Mathieu equation, there are periodic solutions of period either 2*π* or 4*π* on each stability boundary on the parameter plane. Are there periodic solutions (of possibly other periods) for each stability boundary for the asymmetric Mathieu equation? In this section, we investigate this question numerically. In particular, we seek curves on the *δ*–*ϵ* parameter plane where periodic solutions exist; and anticipate that these curves will start from the *δ*-axis at the points from which instability regions emanate (as discussed earlier).

We present here the results obtained for *δ*=0.4201, 0.2689, 0.1867 and 0.1050. The corresponding time periods expected, and obtained, are 4*π*, 10*π*, 6*π* and 8*π*, respectively (note: 10*π* corresponds to *r*=2.5 and not 5 in equation (4.1)). From each of these points on the *δ*-axis, two curves are found to emanate. These curves were obtained using a numerical arc-length based continuation method which is described in appendix A. Results, shown using white lines on the stability diagram, are given in figure 7. The picture, presented here for *ϵ*≥0, is symmetric about the *δ*-axis.

It is clear that, at least for the instability regions considered, stability boundaries correspond to the existence of periodic solutions, even for the unstable regions that are absent for the symmetric Mathieu equation. It seems likely that there are also periodic solutions at all other stability boundaries. In §7, we show theoretically that this is in fact the case.

A few points regarding the numerical search for periodic solutions are now presented.

Since solutions are scaleable, we may assume that the initial conditions at *t*=0 satisfy either *x*_{0}=±1, or , or even . It turns out that for the periodic solutions plotted in figure 7, the solution branch to the left of each instability region has and *x*_{0}=1; while for the solution branch on the right side of each region, *x*_{0}=0 and . The numerical strategy, however, assumes one of these (*x*_{0} or ) to be non-zero (equal to ±1), and lets the numerical routine discover that the other is zero (if in fact it is).

Numerically, we proceed as follows. Starting at *t*=0 with some values of *δ* and *ϵ* (along with, say, fixed and *x*_{0}=1), we numerically integrate forward in time to time *T* (for suitable *T*, which is an even multiple of *π* that we know in advance as indicated earlier), and check that the final state is identical to the initial state. We then iteratively look for a nearby point on the *δ*–*ϵ* plane where these conditions are satisfied again. This well-known procedure, called arc-length based continuation, is described for completeness in appendix A.

Note that these periodic solution branches can (at least potentially) be computed in a small fraction of the time needed to generate, say, figure 5. Thus, as for the symmetric Mathieu equation, finding periodic solutions is an efficient way of computing the stability transition curves for the asymmetric Mathieu equation.

## 6. Theoretical considerations

The asymmetric Mathieu equation (1.1) can be rewritten as(6.1)withandFormally, we observe that equation (6.1) is divergence free, i.e. ∇.** f**=

**0**. This means areas are preserved in the plane.

Let *x*=*R* cos(*ϕ*) and *y*=*R* sin(*ϕ*). Consider a Poincaré map that takes ** x** from the start of a forcing period to the start of the next (i.e. through

*t*=2

*π*). Let (

*R*

_{0},

*ϕ*

_{0}) be an initial point. The Poincaré map sends withfor some as yet unknown continuous functions and . Since solutions to equation (1.1) are scaleable, we can writeSimilarly, scaleability requiresThus, the point (

*R*

_{0},

*ϕ*

_{0}) gets mapped to (

*R*

_{0}

*f*(

*ϕ*

_{0}),

*g*(

*ϕ*

_{0})).

The system is reversible in time. Zero initial conditions lead to zero solutions for all time. These two together imply that non-zero solutions do not become zero within finite time. In turn, this meansWe may, without loss of generality, assume *f*>0.

Consider two nearby initial conditions (*R*_{0}, *ϕ*_{0}) and (*R*_{0}, *ϕ*_{0}+Δ*ϕ*_{0}), as sketched in figure 8. The Poincaré map sends these initial conditions to (*R*_{0}*f*(*ϕ*_{0}), *g*(*ϕ*_{0})) and (*R*_{0}*f*(*ϕ*_{0}+Δ*ϕ*_{0}), *g*(*ϕ*_{0}+Δ*ϕ*_{0})) as sketched in figure 8.

If we think of the entire triangle as composed of initial conditions, then the triangle gets mapped to a triangle because solutions are scaleable (the edges remain straight lines).

Since *f* and *g* are continuous, we writewhere the Δ symbol is now taken to denote ‘small’.

Area preservation now giveswhere we have ignored small quantities of second-order.

We now come to an important point. No matter what finite values we assign to the parameters *δ* and *ϵ*, the functions *f* and *g* depend continuously on them as well as on *ϕ*_{0}. So, if we now change any combination of *δ*, *ϵ* and/or *ϕ*_{0} in any way that we like, *f*(*ϕ*_{0}) always remains finite and non-zero. This means Δ*g*(*ϕ*_{0}) is always non-zero as well, because Δ*ϕ*_{0} is non-zero and positive by choice; it *never* changes sign. It is possible to conclude from a consideration of the case *δ*=*ϵ*=0 and *ϕ*_{0}=*π*/2, i.e.that Δ*g*>0. From here, by continuous changes in parameters and initial conditions, we can arrive at the point of interest to conclude that the absolute value sign may be removed, and sowhich in the limit shows that *g* is differentiable and satisfies(6.2)

It follows that *g*′(*ϕ*_{0})>0 for all *ϕ*_{0}.

We observe from equation (6.2) that if any solution settles down to some stable point *ϕ*^{*}, then *g*(*ϕ*^{*})=*ϕ*^{*} and *g*′(*ϕ*^{*})<1 (the condition for stability of fixed points of iterated scalar maps). This in turn implies that *f*(*ϕ*^{*})>1. Thus, any solution that settles exponentially to some *ϕ*^{*} must grow exponentially in magnitude.

What happens if, instead of a fixed point, *g* has a *k*-cycle, i.e. the *k*th iterate of some *ϕ*^{*} equals itself, or *g*^{k}(*ϕ*^{*})=*ϕ*^{*}? Let *g*(*ϕ*^{*})=*ϕ*_{1}, *g*(*ϕ*_{1})=*ϕ*_{2}, and so on. Then, *ϕ*_{k}=*ϕ*^{*} for a *k*-cycle.

If, in addition, the *k*-cycle is exponentially stable, i.e.(6.3)then it follows thatand the solution grows exponentially in magnitude (the Lyapunov-like exponent is positive).

Consider, now, a gradual change in parameters that causes this unstable point in parameter space to approach a stability boundary. The *k*-cycle (corresponding to *ϕ*^{*}) depends on parameters, and changes gradually as well; it is structurally stable as long as inequality (6.3) holds. Thus, loss of stability can only occur whenat which point we also haveIn other words, an unstable solution (i.e. a growing solution with a positive value of the Lyapunov-like exponent used in this paper; but also a solution where *ϕ*^{*} corresponds to a stable *k*-cycle of the iterated function *g*) can only lose instability (or gain stability; or reach a stability margin) by deforming continuously into a periodic solution as parameters are slowly changed so as to reach a point on a stability boundary.

A stable *k*-cycle in the iterated function *g* implies instability in the system solution. From the corresponding (unstable) point on the parameter plane, moving towards a stable point requires appearance of a periodic solution on the stability boundary.

What we wish to prove, however, is more general. We wish to prove that *every* stability boundary corresponds to the existence of a periodic solution. We will do this by showing that every unstable point in the parameter plane corresponds to the existence of a stable *k*-cycle in the iterated function *g*. Conclusion 6.1 will then be applicable, and the desired result will be established.

Accordingly, we now assume that the solution is unstable, i.e. for some parameter values *δ* and *ϵ*, and referring to equation (3.1),(6.4)This is equivalent to(6.5)where *E* represents expected value, and *n* is sufficiently large that initial transients are not important and the final steady-state behaviour of the iterates of *g*, whether periodic or not, is obtained. From equation (6.2), we havei.e.ThereforeLet(6.6)for some strictly positive number *a*. Moreover, since we assume the system has reached steady-state behaviour, we can also define the variance of ln *g*′ as(6.7)for some *b*>0.

Let *ϕ*_{p}=*g*^{p}(*ϕ*_{0}) for *p*>1. Since *ϕ* is an angle, we can look at its values modulo 2*π*.

The dynamics of the system generates a sequence *ϕ*_{0}, *ϕ*_{1}, *ϕ*_{2}, *ϕ*_{3}, ⋯, which is now bounded (because we look at the values modulo 2*π*). Every bounded sequence has a convergent subsequence (Rudin 1978). Let the subsequence be *ϕ*_{c1}, *ϕ*_{c2}, *ϕ*_{c3}, ⋯, where *c*_{i}>*c*_{j} if *i*>*j*. Let the subsequence converge to *ϕ*^{*}. Then, for any given *ϵ*>0, there exists a finite *M* such that for all *n*>*M*, . We choose two points from the subsequence, not necessarily consecutive, say , with *n*_{2}>*n*_{1}>*M*.

Now consider the function(6.8)Then, applying the chain rule of differentiation,(6.9)Considering *h*′(*ϕ*_{m}) for some *m*>*M* in the subsequence *c*_{1}, *c*_{2}, ⋯, we have

Considering the logarithm of the first-term on the right-hand side, we have (calling it, say, *Z*)For sufficiently larger than (and we are free to choose it so), the central limit theorem applies; in particular, the expected value of *Z* is (see equation (6.6)), and its standard deviation is (see equation (6.7)), which is much smaller. Thus, *Z* is strictly negative, and can be as large as we wish to make it: its exponential is a positive number which can be as small as we like, independent of *ϵ*. It follows thatlies between −1 and 0. In particular, it can be bounded away from 0 by a non-zero amount, such as 1/2, independent of *ϵ*.

Now, consider equation (6.8)Thus, *h*(*ϕ*_{m}) is small; and *h*′(*ϕ*_{m}) is non-zero. By the implicit function theorem, *h* has a zero close to *ϕ*_{m}, say at . This zero corresponds to a stable *k*-cycle of the iterated function *g*.

An unstable solution of the system implies the existence of a stable *k*-cycle of the iterated function *g*.

By Conclusion 6.1, every point on a stability boundary in the parameter plane (or parameter space, if we introduce more parameters) corresponds to the existence of a periodic solution. However, unlike the usual or linear Mathieu (or Hill) equation, the period need not be solely 2*π* or 4*π*, but could be a higher multiple of 2*π*.

## 7. General asymmetric Mathieu equations

More generally, our foregoing results apply to(7.1)where *ϵ*_{1} and *ψ* are additional free parameters (compare with equation (1.1)).

The theoretical results presented above actually hold for the general asymmetric Hill's equation(7.2)with *p*_{1}(*t*+*T*)=*p*_{1}(*t*) and *p*_{2}(*t*+*T*)=*p*_{2}(*t*)∀*t* and for some *T*>0 (*T* can be taken as 2*π* upon scaling time suitably). We are not presently aware of actual physical systems governed by such equations, except for the restricted case we began this paper with.

A few numerically generated stability diagrams for equation (7.1) are now presented. We restrict ourselves to as before along with *ψ*=*π*/4.

A sketch of a subset of the full parameter space is shown in figure 9. Figure 10 shows the stability diagram in the parameter plane *ϵ*_{0}=0, and figure 11 shows the stability diagrams on portions of the two orthogonal planes *ϵ*_{1}=0 and *ϵ*_{0}=0, respectively.

Figure 12 shows the stability diagram on *ϵ*_{1}=0.2 plane and figure 13 shows the stability diagrams on planes *ϵ*_{1}=0.2 and *ϵ*_{0}=0, respectively. Numerical resolution limits the degree of detail that can be trusted in the figures. With more computation time, any of these figures could be regenerated with higher precision and resolution. However, our key point in producing these figures is to emphasize that the white lines, representing numerically obtained periodic solutions, were in fact computed accurate to nine decimal places. Moreover, we did find periodic solutions on *every* stability boundary that we examined, verifying the theoretical results obtained above.

## 8. Conclusions and further work

We have numerically and theoretically studied the asymmetric Mathieu equations, which are strongly nonlinear but conservative, and have scaleable solutions (if *x*(*t*) is a solution, then so is *αx*(*t*) with *α*> 0). We have found that there are infinitely many more instabilities for this system than for the usual Mathieu equation. There are periodic solutions on every stability boundary in the parameter space. The periods of these solutions are not confined to either 2*π* or 4*π*; higher multiples of 2*π* occur. Our theoretical results are also applicable to asymmetric Hill's equations.

Several questions seem interesting and relevant which we have been unable to answer here. For equation (1.1), are there in fact infinitely many instability regions of strictly non-zero width emanating from any finite interval on the *δ*-axis? How do the widths of these regions, say for small *ϵ* and for an *m*:*n* resonance, depend on *m*, *n* and *ϵ*? Given an arbitrary point (*δ*, *ϵ*) with *ϵ*>0, does every open set in the parameter plane that contains this point also contain an unstable point? Given some small damping, how much of which instability regions will survive? We hope that future work may shed some light on these issues.

## Footnotes

- Received July 20, 2005.
- Accepted November 29, 2005.

- © 2006 The Royal Society