0
Research Papers

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

[+] Author and Article Information
Ashu Sharma

Department of Mechanical Engineering,
Auburn University,
Auburn, AL 36849
e-mail: azs0111@auburn.edu

S. C. Sinha

Life Fellow ASME
Alumni Professor Emeritus,
Department of Mechanical Engineering,
Auburn University,
Auburn, AL 36849
e-mail: ssinha@eng.auburn.edu

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

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.

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 [2123] 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.

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

(1)x¨+dx˙+(a+b1cosω1t+b2cosω2t)x=0

where ω1 and ω2 are the two parametric frequencies such that ω1/ω2 is an irrational number (incommensurate); a,b1,b2,andd 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ω1)2/4 and a=(k2ω2)2/4, where k1,k2=0,1,2,. Johnson and Moser [29] showed that in the spectral gaps, the rotation number (α) of an almost periodic function is related to the frequency module (M), 2αM. 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 ab plane originate from a axis at a=(k1ω1±k2ω2)2/4;k1,k2=0,1,2,... where main regions of instability arise from a=ω12/4(k1=1,k2=0) and a=ω22/4(k1=0,k2=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: ω1=1andω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 (ω1=1,ω2=2andω3=3) 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.

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

(2)x˙=A(t)xx(0)=x0

where xn, t+, A(t) is an n×n periodic matrix with the principal period T. Let Φ(t) be the state transition matrix such that it satisfies the above linear equation with Φ(0)=I; then the solution of Eq. (2) can be written asDisplay Formula

(3)x(t)=Φ(t)x00tT

For tT, it can be calculated byDisplay Formula

(4)x(t+sT)=Φ(t)Φs(T)x00tTs=1,2,3,

where Φ(T) is the Floquét transition matrix (FTM) or the monodromy matrix. The stability criteria for periodic systems depend upon the eigenvalues of Φ(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 asDisplay Formula

(5)Φ(t)=L(t)eCtL(t)n×n,Cn×nt 0

orDisplay Formula

(6)Φ(t)=Q(t)eRtQ(t)n×n,Rn×nt 0

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 coefficientsDisplay Formula

(7)z˙=Cz orz˙=Rz,respectively

where C and R are time-invariant matrices.

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

(8)x¨+(a+bcosω1t)x=0

where a=ωn is the natural frequency of the system (with b=0) and ω1 is the parametric excitation frequency. In ab plane, stability boundaries cross b=0 line or originate fromDisplay Formula

(9)a=(Kω12)2K={0,2,4,Tperiodic1,3,5,2Tperiodic

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

(10)x˙=A(θ,λ)x;dθdt=ωx(0,λ)=x0

where x(θ,λ)n, A(θ,λ) is an n×n matrix that is 2π periodic in θ={θ1,θ2,,θm} and is a continuous function of a set of control parameters λ; ω={ω1,ω2,,ωm}m is the frequency vector of A(θ,λ) and t is the time. When A(θ,λ) has finite (m2) incommensurate frequencies, then A(θ,λ) is quasi-periodic.

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

(11)x¨+dx˙+(a+i=1mbicos(ωit))x=0

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

(12)x˙=[01(a+i=1mbicos(ωit))d]x

where x={x1x2}T and ω={ω1,ω2,,ω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 ε>0, there exists a length of time T(ε) that contains a number, Ta(ε) for which |A(t+Ta(ε))A(t)|<ε. For a given ε, Ta(ε) can be determined by truncating the frequency module of A(ωt,λ) defined as |k1ω1±k2ω2±±kmωm|; ki=0,1,2,andk1=k2==km0. Frequency module can be truncated by fixing the upper limit on ki. The minimum frequency, ω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:Display Formula

(13)ωmin=Min(|k1ω1±±kmωm|);ωmin0 andTa=2π/ωmin

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

(14)x˙=[01(a+i=1mbicos(ω¯it))d]x

where ωi in Eq. (12) has been approximated by ω¯i such that ω¯i are integral multiples of ωmin, the minimum frequency of the approximate system. From Eq. (13), this translates into the principal period Ta=2π/ω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;i=1,...,m should be similar to Eq. (9). Thus, the stability boundaries in the approximate system would originate fromDisplay Formula

(15)a=(Kωmin2)2;K={0,2,4,Taperiodic1,3,5,2Taperiodic

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]:Display Formula

(16)Tr(Φ(Ta,λ))=±2(undampedsystem,d=0)

andDisplay Formula

(17)Det(I±Φ(Ta,λ))=0(dampedsystem,d0)

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.

Three cases are studied in this paper. It includes two cases in which A(ωt,λ) in Eq. (12) consists of two frequencies and for the last case A(ωt,λ), 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 ω1=π and ω2=7.0. This set is chosen as a possible representative of systems where ω1 and ω2 are relatively apart. In such cases, finding ωmin could possibly involve a large number of k1andk2 in Eq. (13) as compared to frequency pairs with smaller differences between them. Also, in the corresponding periodic cases with excitation frequencies ω1 and ω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 ω1=1.0 and ω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(ωt,λ), in Eq. (12) contains three frequencies: ω1=1.0, ω2=3, and ω3=11. The difference between the frequencies (ω2ω1=0.732051, ω3ω2=1.58457, and ω3ω1=2.31662) are relatively smaller than the case study 1. A number of instability loops are observed in the parametric space.

Case Study 1: ω1=πandω2=7.0.

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

(18)x˙=[01(a+b1cosπt+b2cos7t)d]x

As a first step, a set of minimum frequencies (ω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 ωmin is defined for a set of maximum values of k1 and k2. For example, as shown in Table 1, ωmin=0.168147 is minimum in the range k1,k2=0to28 and it occurs for k1=20 and k2=9. As k1 and k2 are increased further, a frequency smaller than ωmin=0.168147 is obtained for k1=29 and k2=13. This new ω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 ω1=π and ω2=7 are replaced by ω¯1 and ω¯2, respectively, which are the integral multiples of ωmin. Equation (14) reduces to the formDisplay Formula

(19)x˙=[01(a+b1cosω¯1t+b2cosω¯2t)d]x

Columns 6 and 7 in Table 1 show ω¯1 and ω¯2 for various ωmin. For instance, when ωmin=0.168147 (E. No. 4 in Table 1), ω1/ωmin=18.6836, and ω2/ωmin=41.6302. The nearest integers are 19 and 42 and therefore, ω¯1=19ωmin=3.194793 and ω¯2=42ωmin=7.062174. Since the greatest common divisor (GCD) of 19 and 42 is 1, the ωmin remains as a minimum frequency in the approximate system. In case of ωmin=0.106187 (E. No. 5 in Table 1), ω1/ωmin=29.5855, and ω2/ωmin=65.9214, and thus we can set ω¯1=30ωmin=3.18561 and ω¯2=66ωmin=7.008342. As the GCD of 30 and 66 is 6, the minimum frequency of the approximate system has changed to 6ω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 ω¯1=29ωmin=3.079423 and ω¯2=66ωmin=7.008342.

Once an approximate system has been defined (c.f. Eq. (19)), the STM, Φ(t,λ) can be calculated numerically using Floquét theory. Setting b1=b2=b (for convenience) and using mathematica, Φ(t,a,b,d) of the approximate system was computed and stability charts were plotted in the ab 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(ωt,λ) as the L-P transformation, the change of variable x=P(ωt,λ)z reduces Eq. (10) to the following system [37,38]:Display Formula

(20)z˙=A0(λ)z

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

The Undamped System (d=0).

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 ω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 ω1=π, the stability boundaries of main parametric instability region arise from a=ω12/4=2.46740, while for ω2=7.0, these arise from a=ω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. 13) leads us to the following observations.

For a given ωmin, say ωmin=0.274334, a number of Ta and 2Ta stability boundaries stem from the a -axis in the ab plane, as anticipated (see Fig. 1). Ta and 2Ta unstable regions are either due to primary frequencies (ω¯1 and ω¯2) or due to various combinations of ω¯1 and ω¯2. From Fig. 13, it could be observed that as ω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 ω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 ω1=π and at a=12.25 for ω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ω1±k2ω2)2/4;k=0,1,2, and the main instability regions stem from a=2.46740(k1=1andk2=0) and a=12.25(k1=0andk2=1), which are identical to the Floquét theory. In fact, this expression is same as Eq. (15) if Kωmin is replaced by k1ω1±k2ω2. However, since k1 and k2 are finite by definition of ω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 ω1=π (see Fig. 3). For ω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=2andk2=0.

The selection of ω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=ω12/4 and a=ω22/4, a convergence study of these bifurcation points is used as a guideline in the selection of ωmin. Figure 4 demonstrates the convergence of bifurcation points of the main instability regions as the approximate period, Ta increases (ωmin decreases). The y -axis denotes the difference between the exact (ae) and approximate (aa) bifurcation points, aeaa and the x -axis is represented by the expression, log10(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 ωmin, the accuracy of the proposed method would increase. However, as ωmin is made smaller, Ta increases and hence it would require longer computation time. Keeping the computational time reasonable, ωmin=0.0442270 is used in further investigations. This results in Ta=142.067, ω¯1=71ωmin=3.140117, and ω¯2=158ωmin=6.987866.

Poincaré maps and spectral analysis.

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π/ω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π/ω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 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), 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 ωmin/2 and they represent the influence of other harmonics on the stability boundaries. From columns 6 and 7, it is observed that faafoe, 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 ω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 C1 (a=3.9792 and b=2.5) in the stability chart corresponding to ωmin=0.0619600 (see Fig. 2). Point C1 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 ω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.97649andb=2.5; see Fig. 3), which is close to point C1, 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 C1 are close to foe suggesting that both minimum frequencies, i.e., ωmin=0.0619600and0.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 ab 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ω1±k2ω2)2/4, the frequency content of the boundaries may be expressed as (±k1ω1±k2ω2)/2. With this observation, the general forms of faa and foe listed in columns 6 and 7 (Table 2) may be expressed asDisplay Formula

(21)faafoe(±k1ω1±k2ω2)2(±k1ω¯1±k2ω¯2)2

For instance, point D (a=5.5560andb=4.5) near the stability boundary of R4 (see Fig. 3) has faa=2.27769 and foe=2.28640 (see Fig. 9). Then, faafoe((3ω1+2ω2)/2)=2.2876 suggests that unstable region, R4 is formed due to frequency combination (3ω1+2ω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 ω¯1andω¯2 and their number increases with smaller and smaller ω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 ω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, π 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 π 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.

Generalized Hill’s infinite determinants method.

As discussed earlier, the frequency contents of the stability boundaries can be expressed as (±k1ω1±k2ω2)/2 (c.f. Eq. (21)), implying that the corresponding solution may be written in the form of a generalized Fourier series as Display Formula

(22)x(t)=k1,k2=0+j[Hk1k2cos(k1ω1±k2ω22)t+Lk1k2sin(k1ω1±k2ω22)t]

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(·) and cos(·) 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×41 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ω12ω2)/2, ω1/2, 3ω1/2, and (2ω2ω1)/2, respectively. As the original quasi-periodic system is approximated by a periodic system with excitation frequencies ω¯1=3.140117 and ω¯2=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 ωmin/2. Using Eq. (21), these frequencies may also be expressed as (5ω12ω2)/2, ω1/2, 3ω1/2, and (2ω2ω1)/2, respectively. These are identical to the frequency combinations obtained for the original quasi-periodic system.

The Damped System (d0).

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 ω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=0to15andb=0to4.

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 (ω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.,  through Ĥ (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=2.467b=0.315 and d=0.1) and point B̂ (a=2.467b=0.313 and d=0.1) agree with each other. The solutions (approximate and exact) are unstable and stable at points  and B̂, 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̂ (a=11.230b=2.2 and d=0.1), Ĝ (a=11.250b=2.2 and d=0.1), and Ĥ (a=11.280b=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̂ and Ĥ. However, a disagreement is observed at point Ĝ 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 aaxis. A disagreement at point D̂ (a=4.130b=3.2 and d=0.1) can also be explained similarly. The general behavior of solutions is similar to each other at points Ĉ (a=4.120b=3.2 and d=0.1) and Ê (a=4.140b=3.2 and d=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 ωmin=0.0442270ω¯1=3.140117π) lies almost on the stability boundary of the main instability region of periodic system with excitation frequency ω1=π, 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.

Case Study 2: ω1=1.0 and ω2=(1+5)/2.

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

For this value of ωmin, the approximate system can be defined by replacing ω1=1.0 and ω2=(1+5)/2 by ω¯1=1.00310580 and ω¯2=1.61611490 (c.f. Table 3), and thus it takes the formDisplay Formula

(23)x˙=[01(a+b1cos1.00310580t+b2cos1.61611490t)d]x

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.

The Undamped System (d=0).

Figure 17 shows the stability diagram for the undamped system corresponding to ω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 1.0a5.0 and 0.0b4.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 ab plane and are Ta and 2Ta periodic, respectively. According to Ref. [30], the main instability regions arise from a=0.25(k1=1andk2=0) and a=0.654508(k1=0andk2=1). 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 ω1 (ω2), the main instability region arises from a=0.25 (a=0.654508) and is 2T periodic where T is the period corresponding to excitation frequency ω1 (ω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 ab plane defined by a=0.25to1.0 and b=0.0to1.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 ω1=1.0, ω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 faafoe 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 Damped System (d0).

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=0to5andb=0to1. The approximate boundaries are compared with the exact boundaries computed using the maximal Lyapunov exp