Parametrically excited linear systems with oscillatory coefficients have been generally modeled by Mathieu or Hill equations (periodic coefficients) because their stability and response can be determined by Floquét theory. However, in many cases, the parametric excitation is not periodic but consists of frequencies that are incommensurate, making them quasi-periodic. Unfortunately, there is no complete theory for linear dynamic systems with quasi-periodic coefficients. Motivated by this fact, in this work, an approximate approach has been proposed to determine the stability and response of quasi-periodic systems. It is suggested here that a quasi-periodic system may be replaced by a periodic system with an appropriate large principal period and thus making it suitable for an application of the Floquét theory. Based on this premise, a systematic approach has been developed and applied to three typical quasi-periodic systems. The approximate boundaries in stability charts obtained from the proposed method are very close to the exact boundaries of original quasi-periodic equations computed numerically using maximal Lyapunov exponents. Further, the frequency spectra of solutions generated near approximate and exact boundaries are found to be almost identical ensuring a high degree of accuracy. In addition, state transition matrices (STMs) are also computed symbolically in terms of system parameters using Chebyshev polynomials and Picard iteration method. Stability diagrams based on this approach are found to be in excellent agreement with those obtained from numerical methods. The coefficients of parametric excitation terms are not necessarily small in all cases.

# An Approximate Analysis of Quasi-Periodic Systems Via Floquét Theory OPEN ACCESS

**Ashu Sharma**

**S. C. Sinha**

Contributed by Design Engineering Division of ASME for publication in the JOURNAL OF COMPUTATIONAL AND NONLINEAR DYNAMICS. Manuscript received March 6, 2017; final manuscript received August 17, 2017; published online November 9, 2017. Assoc. Editor: Bogdan I. Epureanu.

*J. Comput. Nonlinear Dynam*13(2), 021008 (Nov 09, 2017) (18 pages) Paper No: CND-17-1104; doi: 10.1115/1.4037797 History: Received March 06, 2017; Revised August 17, 2017

## Abstract

## Introduction

Linear ordinary differential equations with time-periodic coefficients have a wide range of applications in various fields of science and engineering. Some of their applications include structures subjected to periodic loads [1], helicopter rotor blades in forward flight [2], asymmetric rotor bearing systems [3], robots performing repetitive tasks [4], ship dynamics [5], attitude stability of satellites [6], heart rhythms [7] and quantum mechanics [8], among others. These systems are also called “parametrically excited systems” since system parameters in mathematical models are no longer constants but become periodic functions of time. These periodic functions can cause instability known as “parametric resonance.”

In 1883, M.G. Floquét developed a complete theory [9] (now known as the “Floquét theory”) that determines stability and response of linear ordinary differential equations with periodic coefficients. Since then, several numerical, analytical, and symbolic techniques have been developed to compute stability and response of periodic systems using the Floquét theory. Some authors approximated the periodic coefficient matrix (periodic functions) by piecewise constant, linear and quadratic matrix (functions) [10,11] in order to generate approximate solutions. Using Hammond’s improved integration scheme [12], Friedmann et al. [13] employed Hsu’s approach and suggested a numerically efficient method that requires only a single integration pass scheme. Analytical techniques such as perturbation [14] and averaging [15] have been used by a number of authors to determine stability of periodic systems. However, their effectiveness is restricted by the requirement of existence of a small parameter and a generating solution. To overcome the above limitations, Sinha and his associates developed an efficient technique to compute the state transition matrix (STM) of such systems in semi-analytical [16] as well as in symbolic forms [17,18]. Design of control systems to drive arbitrary unstable or chaotic motions to a periodic orbit also utilizes Floquét theory [19,20].

As far as applications are concerned, numerous contributions exist in the literature. In most cases, systems have been modeled by Mathieu or Hill equations because one is able to predict their stability and response by the use of Floquét theory. However, in many cases, these are oversimplified models since in reality the parametric excitation is not periodic but consists of frequencies that are incommensurate, making it quasi-periodic. For example, parametric excitation occurs when a ship is sailing in longitudinal waves, i.e., sea waves are along the length of the ship. Assuming the vertical force due to sea waves as a periodic force, the equation of motion turns out to be a Mathieu/Hill equation. However, sea waves are not periodic in nature; instead, they contain incommensurate frequencies. Like the sea waves, the motion of basilar membrane (BM) in cochlea of the inner ear and the motion of a heart are some other examples where the parametric forcing terms may not be assumed to be periodic. An appropriate modeling of the heart can be helpful in detecting any irregularities/diseases. The dynamics of BM inside cochlea of inner ear has been a mystery for a long time. Later, some authors [21–23] have established the fact that BM undergoes parametric excitation and in case of single tone (frequency) stimuli, dynamic analysis could be performed using the Floquét theory. However, the structure of sound (speech signal) is generally complex and to understand the dynamics under complex stimuli, study of two or three tone dynamics has to be studied first. Robles et al. [24] and Ruggero et al. [25] investigated the response of BM to two tone stimuli and established the presence of distortion products (combination of tones). The origin of quasi-periodic solution can also be traced back to Hopf bifurcation of a fixed point that leads to a periodic solution (limit cycle) whose subsequent bifurcation, generally known as the secondary Hopf bifurcation (or Neimark–Sacker bifurcation), gives birth to either a periodic or a quasi-periodic motion. In order to study the stability of resulting quasi-periodic motion, one has to construct the variational equation about this motion, which gives rise to a quasi-periodic (QP) Hill equation.

The control of instabilities in these systems is as much important as the determination of stable/unstable regions. In order to design a proper control system, it is imperative to take into account the quasi-periodic behavior of these parametrically excited systems. Control systems are also used during heart surgery to reduce the motion of a part of beating heart over which surgery has to be performed [26,27]. All these problems require a theoretical framework to understand their dynamics and control, but unfortunately there is no complete theory for linear differential equations with quasi-periodic coefficients. In the following, an attempt has been made to provide an approximate methodology that can be used for a general class of quasi-periodic systems arising in engineering applications.

## Quasi-Periodic Hill Equation and Related Work

One of the simplest forms of linear differential equations with quasi-periodic coefficients can be represented by a modified damped Mathieu/Hill equation given by

where $\omega 1$ and $\omega 2$ are the two parametric frequencies such that $\omega 1/\omega 2$ is an irrational number (incommensurate); $a,\u2009\u2009b1,\u2009\u2009b2,\u2009and\u2009\u2009d$ are system parameters; and $t$ is the time. Stability of this type of system has been studied by various authors in the past. In 1980, Davis and Rosenblat [28] studied the QP Hill equation using multiple scale technique and computed the stability boundaries. According to them, stability boundaries arise from $a$ axis for two sets of families, $a=(k1\omega 1)2/4$ and $a=(k2\omega 2)2/4$, where $k1,k2=0,1,2,\u2026$. Johnson and Moser [29] showed that in the spectral gaps, the rotation number $(\alpha )$ of an almost periodic function is related to the frequency module $(M)$, $2\alpha \u2208M$. This relationship is valid in the entire parametric space. Since rotation number is independent of the particular solution and is equal to $a$ when $b1=b2=0$, it has been used by several authors to plot the stability diagrams. Zounes and Rand [30] used the results of Johnson and Moser [29] and investigated the QP Hill equation both numerically and analytically. First, the stability diagrams were developed by direct numerical integration and computation of Lyapunov exponents. Later in their paper, expressions for transition curves were derived using regular perturbation and a technique similar to Hill’s method of infinite determinants and a good agreement between analytical and numerical solutions was reported. According to them, regions of instability in $a\u223cb$ plane originate from $a$ axis at $a=(k1\omega 1\xb1k2\omega 2)2/4;\u2003k1,k2=0,1,2,...$ where main regions of instability arise from $a=\omega 12/4\u2003(k1=1,\u2009\u2009\u2009k2=0)$ and $a=\omega 22/4\u2003\u2009(k1=0,\u2009\u2009\u2009k2=1)$. They also mentioned that the higher order resonances are more affected by the damping as compared to the lower order resonances. Broer and Simó [31] numerically explored the QP Hill equation with two parametric frequencies: $\omega 1=1\u2009\u2009and\u2009\u2009\omega 2=(5+1)/2$. They numerically examined the stability boundaries for relatively large values of parameter (as compared to Ref. [30]) using maximal Lyapunov exponent and rotation number. In 2011, Puig and Simó [32] investigated the stability boundaries in a QP Hill’s equation with three parametric frequencies $(\omega 1=1,\u2009\u2009\omega 2=2\u2009\u2009and\u2009\u2009\omega 3=3\u2009)$ and plotted the boundaries using rotation number and maximal Lyapunov exponent.

Perturbation and averaging techniques require existence of small parameters and generating solutions, while the application of the Hill’s infinite determinants method to QP Hill equation becomes progressively tedious with higher order systems and with larger number of terms in the generalized Fourier series. Moreover, it is known that the Hill’s infinite determinants do not converge in all cases [1,33]. Stability charts can be plotted using rotation number and Lyapunov exponents; however, they require extensive computational power and do not provide analytical insight. To overcome the shortcomings of the above techniques and due to the unavailability of a rigorous theory for quasi-periodic systems, an approximate technique has been presented in this paper. The basic premise is based on the suggestion that a quasi-periodic system may be approximately replaced by a periodic system with an appropriate large principal period and thus making it suitable for an application of the Floquét theory. Apart from establishing the proof of concept via numerical integration, a symbolic technique is also presented that is accurate for a wide range of system parameters. First, some basic results of Floquét theory are summarized in the following.

## Some Results From Floquét Theory

Floquét theory predicts the stability and response of linear ordinary differential equations with periodic coefficients [34]. Consider

where $x\u2208\mathbb{R}n$, $t\u2208\mathbb{R}+$, $A(t)$ is an $n\xd7n$ periodic matrix with the principal period $T$. Let $\Phi (t)$ be the state transition matrix such that it satisfies the above linear equation with $\Phi (0)=I$; then the solution of Eq. (2) can be written as

For $t\u2265T$, it can be calculated by

where $\Phi (T)$ is the Floquét transition matrix (FTM) or the monodromy matrix. The stability criteria for periodic systems depend upon the eigenvalues of $\Phi (T)$, called the *Floquét multipliers* and the system is stable if all *Floquét multipliers* lie on or inside the unit circle; otherwise, it is unstable.

According to the Lyapunov–Floquét theorem [35], STM can be expressed as

where $L(t)$ is $T$ periodic and $Q(t)$ is $2T$ periodic. $L(t)$ and $Q(t)$ are known as Lyapunov–Floquét (L-F) transformations. Using the transformations, $x=L(t)z$ or $x=Q(t)z$, Eq. (2) can be reduced to equations with constant coefficients

A specific form of Eq. (2) is the well-known Mathieu equation given by

where $a=\omega n$ is the natural frequency of the system (with $b=0$) and $\omega 1$ is the parametric excitation frequency. In $a\u223cb$ plane, stability boundaries cross $b=0$ line or originate from

## Proposition and a General Methodology

The general linear ordinary differential equations with quasi-periodic coefficients may be represented by

where $x(\theta ,\lambda )\u2208\mathbb{R}n$, $A(\theta ,\lambda )$ is an $n\xd7n$ matrix that is $2\pi $ periodic in $\theta ={\theta 1,\theta 2,\u2026,\theta m}$ and is a continuous function of a set of control parameters $\lambda $; $\omega ={\omega 1,\omega 2,\u2026,\omega m}\u2208\mathbb{R}m$ is the frequency vector of $A(\theta ,\lambda )$ and $t$ is the time. When $A(\theta ,\lambda )$ has finite ($m\u22652$) incommensurate frequencies, then $A(\theta ,\lambda )$ is quasi-periodic.

A simple form of Eq. (10) is the damped QP Hill equation given by

Equation (11) can be rewritten in the state space form as

where $x={x1x2}T$ and $\omega ={\omega 1,\omega 2,\u2026,\omega m}$ is the frequency basis of the coefficient matrix.

The stability of Eq. (11) (or Eq. (12)) cannot be investigated by Floquét theory since the “principal period” for such a system tends to infinity. However, an approximate period can always be defined such that for every $\epsilon >0$, there exists a length of time
$T(\epsilon )$ that contains a number, $Ta(\epsilon )$ for which $|A(t+Ta(\epsilon ))\u2212A(t)|<\epsilon $. For a given $\epsilon $, $Ta(\epsilon )$ can be determined by truncating the frequency module of $A(\omega \u2009t,\lambda )$ defined as
$|k1\omega 1\xb1k2\omega 2\xb1\cdots \xb1km\omega m|$; $ki=0,1,2,\u2026\u2003and\u2003\u2009k1=k2=\cdots =km\u22600$. Frequency module can be truncated by fixing the upper limit on $ki$. The *minimum frequency*, $\omega min$ in the truncated frequency module can be chosen to define the approximate principal period, $Ta$ of the parametric quasi-periodic term. Thus, the following parameters are defined as:

$\omega min$ depends on a particular set of values of $ki$ used in Eq. (13). For example, if the coefficient matrix, $A(\omega \u2009t,\lambda )$ in Eq. (12) consists of two incommensurate frequencies $\omega 1=\pi $ and $\omega 2=7.0$ in the frequency basis, then Eq. (13) can be written as $Min(|k1\omega 1\xb1k2\omega 2|)$. If $k1$ and $k2$ vary from 0 to 14, then $\omega min=0.274334$ and $Ta=22.9034$. This $\omega min$ is observed when $k1=9$ and $k2=4$ but is a minimum up to the maximum value of $k1,k2=19$. Once $\omega min$ and $Ta$ have been calculated, the original quasi-periodic system, Eq. (12), is replaced by the following periodic system:

where $\omega i$ in Eq. (12) has been approximated by $\omega \xafi$ such that $\omega \xafi$ are integral multiples of $\omega min$, the minimum frequency of the approximate system. From Eq. (13), this translates into the principal period $Ta=2\pi /\omega min$. Due to the periodic nature of approximate system, its stability and response can be determined by the Floquét theory (c.f. Sec. 3) and it is expected that the solution of Eq. (14) is the approximate solution of the original quasi-periodic system (12).

As the approximate system is a periodic system with the principal period $Ta$; $Ta$, and $2Ta$ stability boundaries can be obtained like the Mathieu equation as described in Sec. 3 and hence expressions for $a$ when $bi=0;\u2002i=1,...,m$ should be similar to Eq. (9). Thus, the stability boundaries in the approximate system would originate from

Using $Ta$ as the principal period of the system, the STM of the approximate system can be computed and subsequently stability chart can be plotted using the following expressions [36]:

The development of such an approximate theory for quasi-periodic systems also allows one to construct Lyapunov–Perron (L-P) transformation matrices that reduce the linear quasi-periodic systems to systems with time-invariant coefficients. This would be presented in a forthcoming paper.

## Computational Details and Stability Charts for Some Typical Cases

Three cases are studied in this paper. It includes two cases in which $A(\omega \u2009t,\lambda )$ in Eq. (12) consists of two frequencies and for the last case $A(\omega \u2009t,\lambda )$, contains three frequencies. The two frequency cases are the simplest examples of QP Hill equations. In the first case study, the two frequencies are selected as $\omega 1=\pi $ and $\omega 2=7.0$. This set is chosen as a possible representative of systems where $\omega 1$ and $\omega 2$ are relatively apart. In such cases, finding $\omega min$ could possibly involve a large number of $\u2009k1\u2009and\u2009k2$ in Eq. (13) as compared to frequency pairs with smaller differences between them. Also, in the corresponding periodic cases with excitation frequencies $\omega 1$ and $\omega 2$, the main parametric resonance regions would be away from each other and the other unstable regions in between main parametric resonances (due to combination of frequencies) will also be apart making the stability chart possibly less complicated. For the second case study $\omega 1=1.0$ and $\omega 2=(1+5)/2$ (known as golden ratio), where the difference between the two frequencies is 0.618034, much smaller as compared to 3.85841 in case study 1. Existence of an *instability loop* in one of the main instability regions is the typical feature of this case. This case has also been investigated by Broer and Simó [31] using rotation number and maximal Lyapunov exponent. The results obtained from the proposed method have been compared with those given in Ref. [31]. In the third case, the coefficient matrix, $A(\omega \u2009t,\lambda )$, in Eq. (12) contains three frequencies: $\omega 1=1.0$, $\omega 2=3$, and $\omega 3=11$. The difference between the frequencies ($\omega 2\u2212\omega 1=0.732051$, $\omega 3\u2212\omega 2=1.58457$, and $\omega 3\u2212\omega 1=2.31662$) are relatively smaller than the case study 1. A number of instability loops are observed in the parametric space.

The damped QP Hill equation, Eq. (11) reduces to Eq. (1), which can be represented in the state space form using Eq. (12) as

As a first step, a set of minimum frequencies ($\omega min$) are calculated using Eq. (13) over a range of values of $k1$ and $k2$ and are listed with their corresponding approximate periods ($Ta$) in Table 1. It should be observed that $\omega min$ is defined for a set of maximum values of $k1$ and $k2$. For example, as shown in Table 1, $\omega min=0.168147$ is minimum in the range $k1,k2=0\u2009\u2009to\u2009\u200928$ and it occurs for $k1=20$ and $k2=9$. As $k1$ and $k2$ are increased further, a frequency smaller than $\omega min=0.168147$ is obtained for $k1=29$ and $k2=13$. This new $\omega min$ ($=0.106187$) remains minimum up to $k1,k2=48$ as seen in Table 1.

In order to approximate the quasi-periodic system, Eq. (18), by a periodic system, the frequencies $\omega 1=\pi $ and $\omega 2=7$ are replaced by $\omega \xaf1$ and $\omega \xaf2$, respectively, which are the integral multiples of $\omega min$. Equation (14) reduces to the form

Columns 6 and 7 in Table 1 show $\omega \xaf1$ and $\omega \xaf2$ for various $\omega min$. For instance, when $\omega min=0.168147$ (E. No. 4 in Table 1), $\omega 1/\omega min=18.6836$, and $\omega 2/\omega min=41.6302$. The nearest integers are 19 and 42 and therefore, $\omega \xaf1=19\omega min=3.194793$ and $\omega \xaf2=42\omega min=7.062174$. Since the greatest common divisor (GCD) of 19 and 42 is 1, the $\omega min$ remains as a minimum frequency in the approximate system. In case of $\omega min=0.106187$ (E. No. 5 in Table 1), $\omega 1/\omega min=29.5855$, and $\omega 2/\omega min=65.9214$, and thus we can set $\omega \xaf1=30\omega min=3.18561$ and $\omega \xaf2=66\omega min=7.008342$. As the GCD of 30 and 66 is 6, the minimum frequency of the approximate system has changed to $6\omega min$, which implies that the period of the approximate system is not the same as the approximate period of the quasi-periodic system that is computed by truncating the frequency module of coefficient matrix. In order to avoid this situation, we can change one of the integers or both such that the GCD becomes 1. In this particular case, the frequencies of the approximate system could be represented by $\omega \xaf1=29\omega min=3.079423$ and $\omega \xaf2=66\omega min=7.008342$.

Once an approximate system has been defined (c.f. Eq. (19)), the STM, $\Phi (t,\lambda )$ can be calculated numerically using Floquét theory. Setting $b1=b2=b\u2009$ (for convenience) and using mathematica, $\Phi (t,a,b,d)$ of the approximate system was computed and stability charts were plotted in the $a\u223cb$ plane utilizing Eqs. (16) and (17). All computations in this work were performed on a laptop computer with a 2.50 GHz i7-4710MQ 4-core processor and 24 GB of RAM.

*A note of interest*: As indicated in Sec. 3, a linear periodic system can be reduced to a system with constant coefficients by using L-F transformations. A similar result, known as Lyapunov–Perron (L-P) transformations, exists for quasi-periodic systems that reduce the linear quasi-periodic systems to systems with constant coefficients. Denoting $P(\omega t,\lambda )$ as the L-P transformation, the change of variable $x=P(\omega t,\lambda )\u2009z$ reduces Eq. (10) to the following system [37,38]:

where $A0(\lambda )$ is a time-invariant matrix. It should be noted that the frequency basis of L-P transformation, $P(\omega t,\lambda )$, is same as the coefficient matrix, $A(\omega t,\lambda )$. Since in the proposed methodology, the minimum frequency $\omega min$ is computed by truncating the frequency module of $A(\omega t,\lambda )$, the frequency module of the $P(\omega t,\lambda )$ will be same as the truncated frequency module of $A(\omega t,\lambda )$. Column 5 in Table 1 shows the total number of frequencies that could be present in $P(\omega t,\lambda )$ for a particular $\omega min$ and it can be calculated using $((2j+1)m\u22121)/2$ where, $m$ is the number of frequencies in the frequency basis and $j$ is the largest value of $ki;\u2003i=1,2,\u2026$.

In this section, the undamped QP equation is investigated and stability charts are plotted for various values of minimum frequencies in order to ascertain the largest value of $\omega min$ that would yield stability boundaries up to a desired accuracy. As a guideline, it is to be noted that for the periodic case of $\omega 1=\pi $, the stability boundaries of main parametric instability region arise from $a=\omega 12/4=2.46740$, while for $\omega 2=7.0$, these arise from $a=\omega 22/4=12.25$ and both are $2T$ periodic. This is easily verified by Eq. (9). For the QP Hill equation as well, Zounes and Rand [30] reported that stability boundaries of main instability regions arise from same values of $a$. A careful examination of stability charts (see Fig. 1– 3) leads us to the following observations.

For a given $\omega min$, say $\omega min=0.274334$, a number of $Ta$ and $2Ta$ stability boundaries stem from the $a$ -axis in the $a\u223cb$ plane, as anticipated (see Fig. 1). $Ta$ and $2Ta$ unstable regions are either due to primary frequencies ($\omega \xaf1$ and $\omega \xaf2$) or due to various combinations of $\omega \xaf1$ and $\omega \xaf2$. From Fig. 1–3, it could be observed that as $\omega min$ decreases, the number of unstable regions increases in the parametric space and the prominent regions of instability become clearly visible. The bifurcation points on $a$ -axis are $Ta$ and $2Ta$ periodic and can be computed using Eq. (15). It is observed that the two main instability regions stem approximately from $a=2.5$ and $a=12.50$.

For $\omega min=0.0442270$, the main instability regions start from $a=2.46508$ and $a=12.2076$ and are $2Ta$ and $Ta$ periodic, respectively. In periodic systems, the main parametric resonance (2:1) is $2T$ periodic and occurs at $a=2.46740$ for $\omega 1=\pi $ and at $a=12.25$ for $\omega 2=7.0$ (c.f. Eq. (9)), which are close to the values obtained from Eq. (15). Using the concept of winding number, Zounes and Rand [30] reported that all stability boundaries arise from $a=(k1\omega 1\xb1k2\omega 2)2/4;\u2009k=0,1,2,\u2026$ and the main instability regions stem from $a=2.46740$$(k1=1\u2009\u2009and\u2009\u2009k2=0\u2009)$ and $a=12.25$$(k1=0\u2009\u2009and\u2009\u2009k2=1\u2009)$, which are identical to the Floquét theory. In fact, this expression is same as Eq. (15) if $K\omega min$ is replaced by $k1\omega 1\xb1k2\omega 2$. However, since $k1$ and $k2$ are finite by definition of $\omega min$, it is a subset of results given in Ref. [30]. A prominent instability region is also observed near the $T$ periodic unstable region (1:1) of periodic system that emanates from $a=9.86960$ for the excitation frequency $\omega 1=\pi $ (see Fig. 3). For $\omega min=0.0442270$, this instability region is $Ta$ periodic and it stems from $a=9.86033$ that is close to the exact bifurcation point $a=9.86960$ calculated using Zounes and Rand [30] expression with $k1=2\u2009\u2009and\u2009\u2009k2=0$.

The selection of $\omega min$ could involve a number of trials but this problem is circumvented by performing a convergence study of the bifurcation points. Since the main instability regions in the quasi-periodic system arise from $a=\omega 12/4$ and $a=\omega 22/4$, a convergence study of these bifurcation points is used as a guideline in the selection of $\omega min$. Figure 4 demonstrates the convergence of bifurcation points of the main instability regions as the approximate period, $Ta$ increases ($\omega min$ decreases). The $y$ -axis denotes the difference between the exact $(ae)$ and approximate $(aa)$ bifurcation points, $ae\u2212aa$ and the $x$ -axis is represented by the expression, $\u2009log10(Ta/Ta1)$ where $Ta1$ is the smallest $Ta$ in Table 1, i.e., 2.0. A nonuniform convergence is observed. It is expected that with a decrease in $\omega min$, the accuracy of the proposed method would increase. However, as $\omega min$ is made smaller, $Ta$ increases and hence it would require longer computation time. Keeping the computational time reasonable, $\omega min=0.0442270$ is used in further investigations. This results in $Ta=142.067$, $\omega \xaf1=71\omega min=3.140117$, and $\omega \xaf2=158\omega min=6.987866$.

In order to check the accuracy of the stability boundaries, Poincaré maps of approximate solutions (c.f. Eq. (19)) constructed by Floquét theory using $Ta=2\pi /\omega min=142.067$ and numerical (so-called exact) solutions of the original system given by Eq. (18) are compared at some typical points close to the stability boundaries. These points are denoted by A, B, C, D, F, and G in Fig. 3. Approximate and exact (numerical) Poincaré maps constructed at these points are shown in Fig. 5. It is observed that the maps are qualitatively similar; however, there are significant differences between the amplitudes of approximate (left column) and exact (right column) solutions. The differences in amplitudes of approximate and exact solutions are due to slight discrepancies between the location of approximate and exact stability boundaries in the parametric space. For instance, if point B is assumed near the exact boundary by changing parameter $a$ from $a=0.356044$ to $a=0.358719$ keeping $b$ constant ($b=3.6$), a stable exact solution of large amplitude is obtained (see Fig. 6), which is much closer to the approximate Poincaré map for point B (see Fig. 5(ii)). After all, the boundaries constitute the most sensitive portion of the stability chart. The exact ($a,b$) values corresponding to the boundary near these points are determined by computing the maximal Lyapunov exponent via numerical integration. In Table 2, column 3 shows the $b$ values corresponding to these points while columns 4 and 5 show the corresponding $a$ values obtained by the proposed approximate method and by computing the maximal Lyapunov exponents, respectively. Figure 5(vii) shows the Poincaré maps constructed at one of the typical stable points (point H in Fig. 3) and it could be observed that the exact and approximate maps match qualitatively and quantitatively. Maps at other stable points also show similar patterns but are omitted for brevity.

To further ensure that $Ta$ or $2Ta$ stability boundaries constructed using Floquét theory with $Ta=2\pi /\omega min=142.067$ are accurate representations of exact stability boundaries, spectral analysis using discrete Fourier transform (DFT) is performed on approximate and exact solutions. The solutions are generated for points marked A–G in Fig. 3. In Table 2, columns 3 and 4 specify the locations of these points in the parametric space. Column 6 indicates the dominant frequency in the frequency spectrum of approximate solutions ($faa$) near the boundaries that are computed using Floquét theory. As an example, consider point B ($a=0.356044\u2009$ and $b=3.6$) near the main instability region arising from $a=2.46508$. At this point, the approximate dominant frequency $faa=1.56934$. If point B is assumed near the exact stability boundary ($a=0.358719$ and $b=3.6\u2009$), the exact dominant frequency, $foe=1.5708$, which is close to $faa$. Figure 7, shows the frequency spectrum of the solutions at/near point B. Other frequencies in the frequency spectrum are the integral multiples of $\omega min/2$ and they represent the influence of other harmonics on the stability boundaries. From columns 6 and 7, it is observed that $faa\u2248foe$, which implies that the approximate stability boundaries constructed using Floquét theory with $Ta=142.067$ are accurate representations of exact stability boundaries.

As the number of instability regions change with a change in $\omega min$, the characteristic nature of the stability boundaries ($Ta$ periodic or $2Ta$ periodic) representing a particular unstable region in the parametric space can change as well. For instance, consider point C_{1} ($a=3.9792$ and $b=2.5$) in the stability chart corresponding to $\omega min=0.0619600$ (see Fig.
2). Point C_{1} lies near the $Ta$ periodic stability boundary of the prominent instability region that arises from $a=3.68932$. As shown in Figs.
8(i) and 8(ii), the shape of the Poincaré maps constructed at this point for approximate and exact solutions are similar. The DFT of the solutions shows that $faa=1.92076$ (see Fig. 8(iii)). If $\omega min$ is reduced further to 0.0442270, the same unstable region is bounded by $2Ta$ periodic boundaries and it arises from $a=3.70129$. Point C ($a=3.97649\u2009\u2009and\u2003b=2.5$; see Fig. 3), which is close to point C_{1}, is now near the $2Ta$ periodic stability boundary and the approximate dominant frequency, $faa=1.92387$ while near the exact boundary, $foe=1.92855$. The $faa$ of points C and C_{1} are close to $foe$ suggesting that both minimum frequencies, i.e., $\omega min=0.0619600\u2009\u2009and\u2009\u20090.0442270$ could be used in the proposed method.

Discrete Fourier transform is also helpful in identifying the frequency combinations that cause the parametric resonances. Seven prominent instability regions, R1–R7, as shown in Fig. 3, are investigated and the dominant frequencies at/near points A–G are identified and shown in columns 6 and 7 of Table 2. The error depends upon the proximity of a point to the exact boundary and the quality of the DFT software. For periodic systems, the periodicity of the bifurcation points on the $a$ -axis in $a\u223cb$ plane is same as the periodicity of the stability boundaries originating from them. In the case of quasi-periodic systems, since the bifurcation points are given by Ref. [30] $a=(k1\omega 1\xb1k2\omega 2)2/4$, the frequency content of the boundaries may be expressed as $(\xb1k1\omega 1\xb1k2\omega 2)/2$. With this observation, the general forms of $faa$ and $foe$ listed in columns 6 and 7 (Table 2) may be expressed as

For instance, point D ($a=5.5560\u2009\u2009and\u2009\u2009b=4.5\u2009$) near the stability boundary of R4 (see Fig. 3) has $faa=2.27769$ and $foe=2.28640$ (see Fig. 9). Then, $faa\u2248foe\u2248((\u22123\omega 1+2\omega 2)/2)$$=2.2876$ suggests that unstable region, R4 is formed due to frequency combination $(\u22123\omega 1+2\omega 2)$. The parametric resonances due to combination of excitation frequencies as shown in Table 2 are different from combination resonances in periodic systems with multiple degrees-of-freedom where they occur when the excitation frequency is a combination of natural frequencies of different modes of the system.

In between the prominent instability regions, a bunch of narrow strips of unstable regions are observed between R2 and R3 in Fig.
3. These instability regions arise due to various combinations of frequencies $\omega \xaf1\u2009\u2009and\u2009\u2009\omega \xaf2$ and their number increases with smaller and smaller $\omega min$ (c.f. Eq. (15)). The order of resonance $|k1|+|k2|$ of prominent instability regions may be calculated using column 9 of Table 2. From Fig. 3, it is observed that the width of prominent instability regions decreases with an increase in the order of resonance. These narrow strips of unstable regions are of higher order as compared to the prominent regions of instability. The bifurcation points of these narrow strips may not be close to the exact bifurcation points as opposed to the prominent regions of instability that move closer to exact locations and correct widths with a decrease in $\omega min$. It may be conjectured that with further reduction in *minimum frequency*, these narrow strips of unstable regions will also attain appropriate widths and locations in the parametric space. In Sec. 5.1.2, it is shown that in the presence of damping, these instability regions disappear for smaller values of $b$.

In Fig. 10, instability regions of quasi-periodic system are compared with periodic systems with excitation frequencies, $\pi $ and $7.0$. Although, the main instability regions in the quasi-periodic system originate from same bifurcation points on the $a$ -axis as in the periodic systems (c.f. Sec. 5.1.1), the extensions of the boundaries are a bit different from the periodic systems. For smaller values of $b$ ($b<1$), the boundaries of main instability regions are close to each other; however, quasi-periodic boundaries start deviating as $b$ increases. For the quasi-periodic system, in Fig. 3(ii), instability region R6 starts near $a=10.0$, which is similar for periodic case with excitation $\pi $ as shown in Fig. 10(i). The deviation of quasi-periodic boundaries from periodic boundaries could possibly be the influence of other harmonics on a particular instability region. In the following, a technique similar to the Hill’s method of infinite determinants is used to investigate the influence of harmonics on the stability boundaries.

As discussed earlier, the frequency contents of the stability boundaries can be expressed as $(\xb1k1\omega 1\xb1k2\omega 2)/2$ (c.f. Eq. (21)), implying that the corresponding solution may be written in the form of a generalized Fourier series as

Thus, an approach similar to the Hill’s method of infinite determinants may be applied to compute the stability boundaries for quasi-periodic systems. Substituting Eq. (22) in Eq. (18) for the undamped case and equating the coefficients of $sin(\u2009\xb7\u2009)$ and $cos(\u2009\xb7\u2009)$ leads to a set of linear homogenous algebraic equations in terms of $Hk1k2$ and $Lk1k2$. Approximate analytical expressions for stability boundaries can be computed by setting the determinants of coefficient matrices to zero for a finite value of $j$ in Eq. (22). Zounes and Rand [30] used the same approach to plot the stability diagram of a QP Hill equation. However, in their investigation, the amplitude of parametric forcing term was kept rather small ($b=0.1$). In Fig. 11, the stability boundaries of main instability regions obtained from the proposed method using Floquét theory with $Ta=142.067$ are compared with the stability boundaries computed using the Hill’s approach. It is observed that with an increase in the number of terms in the generalized Fourier series (i.e., $j$), the stability boundaries obtained from Hill’s type approach move toward the boundaries obtained from the proposed method.

Although, the Hill’s type approach provides analytical expressions for the stability boundaries, for larger values of $j$, the dimensions of coefficient matrices increase making this approach computationally inefficient. For example, when $j=4$, the dimension of the coefficient matrix is $41\xd741$ and in addition, the convergence of this method cannot always be guaranteed [1,33]. In order to overcome these shortcomings, a symbolic technique for the computation of STM is presented in Sec. 6, and stability diagrams are obtained.

Nevertheless, the Hill’s approach clearly shows the influence of harmonics on the stability boundaries, which in turn explains the presence of a number of peaks in the frequency spectrum of solutions. For example, Fig. 7 shows the frequency spectrum of solutions generated at point B (Fig. 3(i)) near the stability boundary of instability region R2. A number of peaks corresponding to frequencies 0.853538, 1.5708, 4.71239, and 5.42965 are observed in the DFT of the original equation near the exact boundary (see Fig. 7(ii)) and these frequencies can be expressed as $(5\omega 1\u22122\omega 2)/2$, $\omega 1/2$, $3\omega 1/2$, and $(2\omega 2\u2212\omega 1)/2$, respectively. As the original quasi-periodic system is approximated by a periodic system with excitation frequencies $\omega \xaf1=3.140117$ and $\omega \xaf2=6.987866$, the exact peaks in the frequency spectrum are approximately represented by the frequencies 0.863138, 1.56934, 4.71088, and 5.41851 (see Fig. 7(i)). Since the stability boundaries are $2Ta$ periodic, these frequencies are integral multiples of $\omega min/2$. Using Eq. (21), these frequencies may also be expressed as $(5\omega 1\u22122\omega 2)/2$, $\omega 1/2$, $3\omega 1/2$, and $(2\omega 2\u2212\omega 1)/2$, respectively. These are identical to the frequency combinations obtained for the original quasi-periodic system.

In this section, influence of damping ($d$) on the stability diagram of QP Hill equation (Eq. (18)) is studied. Following the procedure described in Sec. 5.1, stability charts are plotted for the damped system using Eq. (17). Figure 12 shows the stability chart for $\omega min=0.0442270$ and $d=0.1$. It is observed that only some of the prominent instability regions persist in the parametric region defined by $a=0\u2009\u2009to\u2009\u200915$$and\u2009\u2009b=0\u2009\u2009to\u2009\u20094$.

Stability boundaries from the proposed method are compared with the exact boundaries obtained using Lyapunov exponents and Poincaré maps. Since maximal Lyapunov exponent is negative for a stable damped system and positive for an unstable system, the exact stability boundaries can be plotted by detecting the critical values of system parameters where the maximal Lyapunov exponent changes from positive to negative or vice versa. For the unstable region R2, the exact stability boundaries obtained from maximal Lyapunov exponent lie on the approximate boundaries computed using Floquét theory with $Ta=142.067$ ($\omega min=0.0442270$). A slight variation between approximate and exact boundaries can be observed in case of other prominent instability regions. Poincaré maps at some typical points viz., $A\u0302$ through $H\u0302$ (Fig. 12) are constructed and shown in Figs. 13(i)–13(viii). For the unstable region R2, the Poincaré maps of approximate and exact solutions at point $A\u0302$ ($a=2.467$$b=0.315$ and $d=0.1$) and point $B\u0302$ ($a=2.467\u2009\u2009b=0.313$ and $d=0.1$) agree with each other. The solutions (approximate and exact) are unstable and stable at points $A\u0302$ and $B\u0302$, respectively. It is important to note that there is a very small difference between these points and the stability boundary lies in between them. Now, consider Poincaré maps constructed at points $F\u0302$ ($a=11.230\u2009\u2009b=2.2$ and $d=0.1$), $G\u0302$ ($a=11.250\u2009\u2009b=2.2$ and $d=0.1$), and $H\u0302$ ($a=11.280$$b=2.2$ and $d=0.1$). Although the growth/decay rate of solutions is not the same, the qualitative behaviors of exact and approximate solutions match at points $F\u0302$ and $H\u0302$. However, a disagreement is observed at point $G\u0302$ since the exact solution is found to be stable while the approximate solution is unstable. This discrepancy is expected as the approximate bifurcation point, $aa=12.2076$ lies before the exact point, $ae=12.25$ on the $a\u2212axis$. A disagreement at point $D\u0302$ ($a=4.130\u2009\u2009b=3.2$ and $d=0.1$) can also be explained similarly. The general behavior of solutions is similar to each other at points $C\u0302$ ($a=4.120$$b=3.2$ and $\u2009d=0.1$) and $E\u0302$ ($a=4.140\u2009\u2009b=3.2$ and $\u2009d=0.1$).

In Fig. 14, stability boundaries of main instability regions are plotted using three approaches: maximal Lyapunov exponent, proposed method using Floquét theory with $Ta=142.067$ and Hill’s type approach with $j=1$ in Eq. (22). A good agreement among all three approaches is observed near $bcr$, the minimum amplitude of parametric forcing term required for instability. In Hill’s type approach, a better approximation can be achieved by including more terms in the generalized Fourier series (c.f. Eq. (22)). However, no such attempts were made.

Thus far, $b1$ and $b2$ in Eq. (18) have been assumed to have the same values for the sake of convenience. Figure 15 shows the effect of $b2$ on the main instability region, R2 (see Fig. 3). When $b2=0$, the approximate stability boundary constructed using Floquét theory with $Ta=142.067$ (or $\omega min=0.0442270\u21d2$$\omega \xaf1=3.140117\u2248\pi $) lies almost on the stability boundary of the main instability region of periodic system with excitation frequency $\omega 1=\pi $, as anticipated. As the value of $b2$ increases, the instability region R2 moves toward the vertical axis, $a=0$ implying that with the increase in $b2$, more unstable regions are squeezed into the parametric space. It should be noted that by changing the value of $b2$, an unstable point can be stabilized and the bifurcation point can also be moved.

In this case, the frequencies form the “golden ratio” that is well known for its interesting history and irrationality. In contrast to the previous case study, the frequencies $\omega 1$ and $\omega 2$ are smaller in magnitudes and also the difference between them is relatively small, i.e., 0.618034 instead of 3.85841 in case study 1. Due to the proximity of the frequencies, it is anticipated that more instability regions will be squeezed into a smaller region, and thus the stability chart could be quite intricate.

The damped QP Hill equation is similar to Eq. (18) with excitation frequencies 1.0 and $(1+5)/2$. Following the procedure discussed in Sec. 5.1, $\omega min$ is computed over a range of $k1$ and $k2$ and are tabulated in Table 3 along with their corresponding periods, $Ta$. It is important to note that $\omega min=0.0212862$ (E. No. 8 in Table 3) is almost one half of the $\omega min=0.0442270$ (minimum frequency used in case study 1, see E. No. 7 in Table 1) but it is obtained for relatively smaller value of $k1$ and $k2$ ($k1=34$ and $k2=21$). This occurs due to the smaller difference between the frequencies.

Once $\omega min$ has been fixed, the approximate system can be defined using columns 6 and 7 of Table 3. The selection of $\omega min$ was based on a convergence study of the bifurcation points of the main instability regions stemming from $a$ -axis corresponding to $\omega 1=1.0$ and $\omega 2=(1+5)/2$. Figure 16 shows the convergence study of two bifurcation points and it can be observed that for $\omega min<0.0901699$, the difference between the exact ($ae$) and approximate ($aa$) bifurcation points of main instability regions (i.e., $ae\u2212aa$) is relatively small. For $\omega min=0.0557281$, $ae\u2212aa$ is −0.001555 and +0.001552 for regions R1 and R2 (see Fig. 17), respectively. The difference $ae\u2212aa$ decreases with a decrease in $\omega min$; however, this causes $Ta$ to increase and as a compromise, $\omega min=0.0557281$ ($Ta=112.747$) was selected as the minimum frequency for further investigations (see Table 3).

For this value of $\omega min$, the approximate system can be defined by replacing $\omega 1=1.0$ and $\omega 2=(1+5)/2$ by $\omega \xaf1=1.00310580$ and $\omega \xaf2=1.61611490$ (c.f. Table 3), and thus it takes the form

In the following, the STM is computed numerically using Floquét theory with $Ta=112.747$ for undamped and damped systems and stability charts are plotted. For the sake of convenience, all computations are performed by setting $b1=b2=b$.

Figure 17 shows the stability diagram for the undamped system corresponding to $\omega min=0.0557281$ (or $Ta=112.747$). As anticipated, the stability chart is quite complicated as compared to the case study 1. The investigation is restricted to the parametric space of $\u22121.0\u2264a\u22645.0$ and $0.0\u2264b\u22644.0$ and prominent instability regions are marked as R1–R11.

The main instability regions R1 and R2 stem from $a=0.251555$ and $a=0.652957$ in the $a\u223cb$ plane and are $Ta$ and $2Ta$ periodic, respectively. According to Ref. [30], the main instability regions arise from $a=\u20090.25$$(k1=1\u2009\u2009and\u2009\u2009k2=0\u2009)$ and $a=\u20090.654508$$(k1=0\u2009\u2009and\u2009\u2009k2=1\u2009)$. It is to be noted that the approximate bifurcation points calculated using Floquét theory with $Ta=112.747$ are extremely close to the exact bifurcations that are computed using the concept of winding number. In case of periodic system with excitation frequency $\omega 1$ ($\omega 2$), the main instability region arises from $a=\u20090.25$ ($a=\u20090.654508$) and is $2T$ periodic where $T$ is the period corresponding to excitation frequency $\omega 1$ ($\omega 2$). The approximate and exact bifurcation points of other prominent instability regions, marked as R3–R11 in Fig. 17, are listed in the last two rows of Table 4 and a good agreement between them can be observed.

This case has also been studied by Broer and Simó [31] using rotation number and maximal Lyapunov exponent. They numerically plotted the stability diagram in the $a\u223cb$ plane defined by $a=0.25\u2009\u2009to\u2009\u20091.0$ and $b=0.0\u2009\u2009to\u2009\u20091.0$. In Fig. 18, the approximate stability boundaries computed using the Floquét theory with $Ta=112.747$ are compared with the numerically computed boundaries reported in Ref. [31]. The comparison is made at some typical points (marked as black dots) and it is clear that there is an excellent agreement between the two results in the main instability regions.

One of the special features of this case is the presence of an *instability loop* in the main instability region, R2 (see Fig. 17). The *instability loops* or *instability pockets* [39,40] are known to appear in periodic cases where the parametric forcing term has two or more commensurate frequencies. These pockets also exist in Meissner’s equation [41]. Broer and Simó [31] were unable to capture this feature since their investigation was restricted to a smaller parametric region. These loops also exist in other unstable regions (for instance, R6 in Fig. 17).

The accuracies of the stability boundaries were also checked by constructing Poincaré maps of the approximate and exact (numerical solution of Eq. (12) with $\omega 1=1.0$, $\omega 2=(1+5)/2$ and $d=0$) solutions at some typical points (marked A–K in Fig. 17) in the stability chart. The exact and approximate maps were found to be qualitatively similar at all these points. For the sake of brevity, these maps are not shown here.

In order to check the frequency contents of approximate and exact boundaries, DFT analyses were performed at points A–K. The results are listed in Table 4, and it can be observed that $faa\u2248foe$ at all points implying that the boundaries computed using the Floquét theory are accurate representations of the exact boundaries. The dominant frequencies in approximate solutions are further used to identify the frequency combinations causing parametric resonances and the results are also listed in Table 4.

The approximate stability chart obtained from Eq. (23) with $d=0.1$ is shown in Fig. 19. Only a few prominent instability regions remain in the parametric space defined by $a=0\u2009\u2009to\u2009\u20095\u2009\u2009and\u2009\u2009b=0\u2009\u2009to\u2009\u20091$. The approximate boundaries are compared with the exact boundaries computed using the maximal Lyapunov exp