SOLUTION AND STABILITY OF A SET OF PTH ORDER LINEAR DIFFERENTIAL EaUA TIONS WITH PERIODIC COEFFICIENTS VIA CHEBYSHEV POLYNOMIAlS

Chebyshev polynomials are utilized to obtain solutions of a set of pth order linear differential equations with periodic coefficients. For this purpose, the operational matrix of differentiation associated with the shifted Chebyshev polynomials of the first kind is derived. Utilizing the properties of this matrix, the solution of a system of differential equations can be found by solving a set of linear algebraic equations without constructing the equivalent integral equations. The Aoquet Transition Matrix (FTM) can then be computed and its eigenvalues (Aoquet multipliers) subsequently analyzed for stability. Two straightforward methods, the 'differential state space formulation' and the 'differential direct formulation', are presented and the results are compared with those obtained from other available techniques. The well-known Mathieu equation and a higher order system are used as illustrative examples.


INTRODUCTION
The study of systems governed by a-set of ordinary differential equations with periodic coefficients is of great importance in diverse branches of science and engineering.The stability and response under various excitations are the key issues discussed in the vast amount of classical literature available on this subject.Numerous practical applications can be found in the areas of quantum mechanics, dynamic stability of stuctures, circuit theory, systems and control, and dynamics of rotating systems, among others.
In the past, several methods have been used to study the stability of systems with periodic coefficients.These include Hill's method [1,2], the perturbation method [3], the averaging approach [4], and Floquet theory with numerical integration [5].It is well known that both Hill's method and Floquet theory with numerical integration can only be used to determine the stability boundaries as they do not yield closed form solutions for all time.In addition, the former is not computationally convenient for large order systems.The perturbation and averaging methods have their own limitations due to the fact that they can only be applied to systems where the periodic coefficients can be expressed in terms of a small parameter.
A number of authors [6,7,8] have tried to determine the stability and response from an approximate system of equations, which are usually obtained by replacing the elements of the periodic coefficients matrix by piecewise constants or linear functions.In practical application, one can approximate the periodic matrix at most by a series of step functions and compute the transition matrix during one period, which then yields the stability conditions, etc.Such a technique was developed by Hsu [6] and employed by Friedmann et al. [9] for a numerical evaluation of the transition matrix.Although the approach is straightforward, it is only a second-order algorithm at the most.For more accurate solutions, it is necessary to apply higher order numerical schemes such as the Runge-Kutta-Gill method or similar algorithms, and utilize Floquet theory to establish stability conditions.This approach has been adopted by several authors [9,10] in a variety of stability and response problems, and it has been shown by Gaonkar et al. [10] that the use of Hamming's fourth-order predictor-corrector method in a single-pass scheme is very likely the most economical approach.
Recently several studies (e.g.[11][12][13][14]) have been reported where the solutions of both constant-coefficient and time-varying systems are expressed in terms of Chebyshev polynomials.However, the first applications of orthogonal polynomials to differential equations with periodic coefficients were reported by Sinha and Chou [8] and Sinha et al. [7].These applications were limited to second-order scalar equations only.In a later study by Sinha and Wu [15], a general scheme for the stability analysis of a system of second-order equations was presented.The approach was based on the idea that the state vector and the periodic matrix of the system can be expanded in terms of Chebyshev polynomials over the principal period.Such an expansion reduces the original problem to a set of linear algebraic equations, from which the solution in the interval of one period can be obtained.Furthermore, the technique was combined with Floquet theory to yield the transition matrix at the end of one period and provide the stability conditions via an eigen-analysis procedure in addition to providing the solution for all time.Two formulations, one applicable to a set of equations written in state space form and the other suitable for direct application to a set of second-order equations, were presented.In both formulations, the original system of differential equations was converted to a set of integral equations.An error analysis concluded that the suggested schemes not only provide accurate results with rapid convergence, but are also computationally very efficient.In particular, the direct formulation was found to be several times faster than the standard Runge-Kutta type codes.
This procedure was further developed and applied to a large-scale rotorcraft problem [16,17], a boundary value problem in the stability analysis of slender rods [18], and to the optimal control problem of mechanical systems subjected to periodic loadings [19].In particular, it was found by Sinha and Joseph [20] and Sinha et al. [21] that the well-known Liapunov-Floquet (L-F) transformation may be employed which converts the original time-periodic system into an equivalent time-invariant one suitable for the application of standard time-invariant methods of control theory.More recently, this approach has been extended by Sinha and Pandiyan [22], Sinha et al. [23], and Pandiyan et al. [24] to nonhomogeneous and quasilinear systems, in which the time-periodic linear part becomes completely time-invariant through the application of the L-F transformation and the well-known methods of normal forms or averaging can be applied after the transformation to yield more accurate results.
All these developments in Chebyshev polynomial approach, as described above, were basically extensions of the technique originally suggested by Sinha and Wu [15].The differential equations were first rewritten in their equivalent integral forms and the solutions were then obtained by utilizing properties of the operational matrix of integration associated with the shifted Chebyshev polynomials of the first kind.Further, the 'direct approach' of Sinha and Wu [15] is limited to a set of second order equations only.Although the 'direct approach' is much more efficient as compared with the 'state-space method', it involves lengthy manipulations through integration by parts.Further, the solution of sets of larger order equations requires more integrations by parts, and the formula for the algebraic equations quickly becomes unreasonably long.In this paper, a new method is presented which, for the first time, utilizes Chebyshev polynomials to solve a system of pth order linear differential equations with periodic coefficients,where p is an arbitrary positive integer.For this purpose, first the operational matrix of differentiation associated with the shifted Chebyshev polynomials of the first kind is derived.Utilizing the properties of this matrix, it is possible to obtain solutions of a set of pth order equations without constructing the equivalent integral equations.This 'differential formulation' is much more straightforward and efficient when the order of equations is large.Although the initial conditions do not enter as an integration constant as in the 'integral formulation', they may be incorporated into the algebraic equations by taking into account the properties of the differentiation operational matrix.Two methods are presented: one applicable to a set of pth order'equations rewritten in state space (first order) form, and the other directly applicable to the orginal set of pth order equations.It is shown that fewer algebraic equations must be solved in the 'direct method', thus making this approach much more efficient when either n (the number of equations) or m (the number of Chebyshev terms used) is large.The Mathieu equation and a higher order system of two coupled Mathieu equations (which arises when two coupled pendulums have independent vertically oscillating supports) are used as illustrative examples.

Properties of Shifted Chebyshev Polynomials
The Chebyshev polynomials of the first kind are defined by the expression [25,26] and are orthogonal over the interval [-1,1] with respect to the weight function ( These polynomials can also be obtained by the formula T,(t) = cos (76); t = cos(6), -1 =5t =51, r = 0, 1,2, 3. . . s.c.SINHA The shifted Chebyshev polynomials of the first kind are defined in terms of the Chebyshev polynomials of the first kind by using the change of variable, t* = (t + 1)/2, 0 s; t* s; 1.
(4) (5) All properties of T;(t) can be deduced from those of Tr(2t -1).By identifying the first few successive terms, it is also possible to express small powers of t in terms of the shifted Chebyshev polynomials of the first kind as .
The recurrence relation of the shifted Chebyshev polynomials of the first kind can be generated from equation (5) as T;+I(t) = 2(2tl)T;(t) -T;-I(t) (7) while the orthogonality relationships are given by { 0, r *-k I.. 'IT fo T ,(t)T ,(t)w(tjdt ~.~: : : : : : where is the weight function of the shifted Chebyshev polynomials of the first kind.Generally, an arbitrary continuous time function f( t) can be expanded into a shifted Chebyshev series over the interval [0,1] as [27] fit) = L arT;(t); 0 s; t s; 1. r=O (10) The Chebyshev coefficients ar can be obtained from ar ="8 0 w(-r)f(T)T r(T)dT; r = 0, 1,2, 3. . .(11) where w(-r) is the appropriate weight function and This may be written in vector fonn as d .
From the differentiationpropertyof Chebyshevvectors,it can also be shownthat where ( )Tdenotes the transpose of the quantity ( ).The differentiation operational matrix of shifted Chebyshev polynomials of the first kind can also be obtained in an arbitrary interval [t],t2]as shown in Appendix A.2. Since the process of differentiating k times reduces the order of the polynomials by k, the D matrix is a nilpotent matrix of rank m -1 containing zeros on its main diagonal and As discussed in references [15,25,26], the process of integrating a function which has been expanded in m terms of shifted Chebyshev polynomials using the integration operational matrix G (see Appendix B) results in a type offorward difference recurrence procedure in which a loss of information occurs because the resulting expression, instead of having m + 1 terms, is truncated after m terms in order to keep the polynomial vector the same length.In that procedure, the initial conditions enter when the integration constant is added, however, thus keeping the number of Chebyshev terms (and hence the accuracy of the expansion) constant.On the other hand, due to the nilpotent nature of D and DT, the process of differentiating a function which has been expanded in shifted Chebyshev polynomials using the differentiation operational matrix results in a type of backward difference recurrence procedure in which a loss of accuracy occurs because the constant coefficient (which is usually a dominant term in the Chebyshev expansion) is lost and the number of Chebyshev terms decreases by one for each differentiation.Hence, the initial conditions do not enter automatically, and must be added in such a way as to keep the resulting number of algebraic equations the same.
This property of the differentiation operational matrix can be verified easily by successively integrating and differentiating a shifted Chebeshev polynomial vector of length m.Such an operation without the use of operational matrices simply yields the original vector back.However, this is not the case when the integration and differentiation operational matrices G and D are used.As an illustration, consider the expressions 0 0 0 ... 0 It is to be noted that the product of the differentiation (D) and integration (G) operational matrices (respectively the transpose of this product) results in an (m -1) X (m -1) identity matrix along with a zero row (respectively column) and extraneous column (respectively row) rather than an m X m identity matrix.It can also be concluded that integrating and differentiating a polynomial vector (of length m) k times with the use of operational matrices yields a coefficient matrix which contains an (mk) X (m -k) identity matrix along with k zero rows (respectively columns) and extraneous columns (respectively rows).

METHODOF ANALYSIS
Consider an n dimensional linear system of pth order differential equations of the form Lit)y(t) = 0; = ; 0, ...
wherey(t) is an n X 1 position vector (YI(t)Y2(t)...Yn(t))Tandl/t) is the linear differential time-periodic vector operator defined as Wk and Ak(t), k = 0,1,.. .,pare n X n constant and time-periodic matrices, respectively.Each of the time periodic matrices Ak(t) can be written in the most general form as where the functions fi(t) = fi(t + 13 i), I = 1,.. .,rkare periodic with period 13 i and the n X n constant matrices A i, I = 1,.. .,rkcontain the coefficients of these periodic functions.Since the frequencies are commensurate, there exists a positive number T such that q ~13 ~= T for positive integers q~.Note that the minimum value of T for which this relation holds is known as the 'principal period' and is precisely the period of the entire operator trY) = tit + 1).The transformation t = TT normalizes the period of the functions fk(T) to l/q ~and the resulting equation is The new operator with new matrices defined by Wk = k Wk and AiT) = Ak(T + 1) = kAiT), T T k = 0, 1,. ..,p.In the following, the solution of equation ( 20) is obtained in terms of shifted Chebyshev polynomials of the first kind by two different methods via the differentiation operational matrix.
At this stage the lIq i-periodic functionsf~( 7), K = 0, I, A = I,.. "PKare expanded in m terms of shifted Chebyshev polynomials of the first kind with known coefficients d ~Aand the state vector X( 7) is similarly expanded but with unknown coefficients b{.These take the forms m-' f~ (7)  where T; (7) (0 :5 7 :5 1) are the shifted Chebyshev polynomials of the firstkind andXi 7) are the elements of the state vector X (7).Substituting equation (24) into equation ( 22) yields T p;, (7)[WI + QdD pJJpn + T pn (7) ( As a result of the zero row of DT discussed in section 2.2, jj ;n contains pn zero rows (from the Kronecker product operation) which are located in the imth rows where i = 1,..., pn.The corresponding rows of Zpn and 0 may therefore be removed and replaced by the pn initial conditions (cf equation ( 22)) expressed in terms of the Chebyshev coefficients as where ( 29) This guarantees that the initial conditions are satisfied and equation ( 22) is transformed to a set of pnm nonhomogeneous linear algebraic equations for the elements of Bpn given by where Zpn is the modified Zpn matrix appearing in equation ( 28) and go = (00 ...OX~OO ...OX~.. .00 ...OX~n)T contains the elements of the initial condition vector XO.The above equation yields the Chebyshev coefficients similar to the backward recurrence procedure given in reference [25].Finally, the state vector may be determined from AT Since the application of the differentiation operational matrix results in the removal of one of the Chebyshev polynomials, it is, therefore, necessary to increase the number of polynomials used in this process by one in order to provide an equivalent accuracy with the integral formulation [15].It is to be noted that in this formulation, matrix inversion has been completely avoided as far as converting the problem to a state space form is concerned.
It is observed from equation ( 30) that in the 'differential state space formulation', a set of pnm equations must be solved.In the following section.analternate formulation is presented such that the total number of algebraic equations is considerably reduced.

Differential Direct Formulation
Returning to equation (2Q.),each ~the peri~c matrices appe~g in equation ( 21 where YiT) are'the elements of Y(T).Substituting equation (33) in equation ( 20 The initial conditions are thus guaranteed to be satisfied and equation ( 20) is transformed to a set of nm nonhomogeneous linear algebraic equations for the elements of Bn given by where Zn is the modified Zn matrix appearing in equation ( 35 contains the elements of the initial condition vectorsyO,;o,...,y(P-l)o.Equation (37) may be solved for all of the Chebyshev coefficients Bn and the position vector and its derivatives may be determined from or, since i'Jn(T) consists of p blocks of i'~(T) on its main diagonal, the complete state vector X(T) = (y T(T); T(T)' .-y,(p-l)T(T)) T may be obtained from equation (32) where Since the application of the pth power of the differentiation operational matrix results in the loss of p Chebyshev polynomials, it is therefore necessary to increase the number of polynomials used in this process by p in order to provide the same accuracy as in the integral formulation, in contrast to the 'differential state space formulation' in which it is necessary to increase the number of polynomials only by one.However, this approach is much more efficient since the 'differential state space formulation' requires a solution of pnm algebraic equations in contrast to only nm equations in the 'differential direct formulation'.To seethis,the substitutionm ~m + 1is madetopnm andthe substitution m ~m + p is made to nm.Sincethese have the effectof addingpn equationsin both cases, the difference remains the same, i.e. (p-l)nm more equations must be solved in the 'differential state space formulation' than in the 'differential direct formulation.'

Computation of the Floquet Transition Matrix
Floquet theory is the most significant tool in the analysis of stability and response for linear differential equations with periodic coefficients.In this section we briefly summarize some main results and show how the Floquet transition matrix (FfM) can be computed using the aforementioned methods.
Consider a set of linear ordinary differential equations of the form where x( t) is an n vector and A( t) is an n X n periodic matrix having a principle period T, such that A(t + T) = A(t).Based on the Floquet theorem [5], the fundamental matrix 4>(t)can be represented as 4>(t)= P(t)exp(Zt), where P(t) = P(t + T) is a periodic matrix, and Z is a constant matrix.In general, «I»(t)is a solution of equation ( 40 Moreover, for times greater than one period, the solution is obtained as x(t + kT) = 4>(t

+ k1)x(O) = «I»(t)[«I»(1)]"x(O).
The knowledge of F = «1»(1) is significant in the study of periodic differential equations since it provides necessary and sufficient conditions for the stability of linear periodically time-varying systems.The eigenvalues ~kof «1»(1) are called the Floquet multipliers and, in general, they are complex (~k = ~kR + i~kI)' We can also define the characteristic exponents as (Xk:t illk,where (Xk= (lI1)lnl~kl and Ilk= (lI1)tan-I(~kI/~kR)' The stability criteria are related to the characteristic exponents and Floquet multipliers and may be summarized as in the following: (i) The knowledge of «1»(1) at the end of one period is sufficient to predict the stability of the system.The stability criteria for the system can be stated in terms the real part of characteristic exponent (Xkor in terms of the modulus of the Floquet multiplier.A system is stable if all (Xk< 0 or if alll~kl < 1 for k = 1,2...,n.If at least one (Xk> 0 (or at least one I~kl> 1), it is unstable; (ii) Once the fundamental matrix «I»(t)is known in the interval [0,1], the solution for all t > T can simply b~obtained ftom the semigroup property as stated above.The details can be found in reference [28].

([D~F-1Bn) T)T.
The solution for the state vector in both the 'differential state space formulation' and the 'differential direct formulation' is given by equation (32).This is valid in the interval t e [0,1] or T e [0,1].As pointed out previously, the solution can be easily extended for t > T (T > 1) by utilizing the formula

MATHIEU'SEQUATION
We first consider the well-known problem of the Mathieu equation, Q is the product operational matrix associated with COS(2'TTT), and DT, I, and Q are all of dimension m X m due to m number of polynomials used in the Chebyshev expansion.The zero mth and 2mth rows of jj Icorrespond with the rows of Z2 and 0 that are replaced with the 2 initial conditions as in equation ( 29).The resulting system of 2m algebraic equations is of the form of equation ( 30) and may be solved for the Chebyshev coefficients B2 to obtain the solution.
The 'differential direct formulation' (cf section 3.2) reduces to equation (34) with 4'TT2 where all matrices are defined as in equation ( 47) and we have again multiplied through 4'TT2 by 02.The initial conditions (equation (36» replace the (m -l)th and mth rows of Zl and 0 and the resulting system of m algebraic equations in B 1is of the form of equation ( 37), which can also be solved to obtain the solution.
'qle numerical results for the case of a = 1,b = -0.32,and 0 = 3 werecomputedvia both the 'differential state space' and 'differential direct' methods over one period and compared with the results obtained previously by Sinha and Wu [15] for the 'integral formulation' and a Runge-Kutta (DVERK) routine in the IMSL library.For expansions of 10 and 12 terms, the results using the 'differential state space formulation' agree precisely up to 6 significant figures with those of the 'integral formulation' for 10 and 12 term expansions, respectively, and with the Runge-Kutta routine up to five significant figures.As expected, the 'differential direct fonnulation' needed 11 and 13 terms, respectively, for an equivalent accuracy.The familiar stability diagram is presented in Figure 1, where 0 = 2'TT and 67,200 different points in the two-dimensional (a,b) parameter space have been analyzed for stability by computing the eigenvalues of the FfM found by equation ( 42).
Here, the dark regions represent unstable zones (when one eigenvalue is outside the unit circle in the complex plane) and the white regions the neutrally stable ones (when both eigenvalues lie on the unit circle).Since there is no damping present, all of the unstable 'Arnold tongues' [29] intersect the b = 0 line at a = (q'TT)2, q = 1,2,...in this Hamiltonian system; however, the finite resolution does not show this for the one farthest right.It is to be observed that the boundaries between the zones exactly coincide with those reported in [30] and [31] where the b = 0 intersections are at l/4, q = 1,2,...(for 0 = 1) and that the route of destabilization is specific for each tongue, alternating between tangent or saddle-node (when the eigenvalues of the FfM meet and split at + 1) and period doubling (when they meet and split at -1).Because the FfM for a Hamiltonian system is a s.c.SINHA symplectic mapping whose eigenvalues must lie symmetrically opposite both the real axis and the unit circle in the complex plane, these are the only possible routes of destabilization (bifurcations in the corresponding nonlinear equation) for a single degree-of-freedom system with two eigenvalues [32].

A Higher Order System
We next consider the problem of two coupled Mathieu equations of the form each of which has period T = n and parameters ai' bi, and Ci'i = 1,2.This is the model for the linearized system of two coupled pendulums with vertically oscillating supports where ai = g/li + klmi, bi = Ai 0.2/1li'and Ci = -kli,~jm;li' i = 1,2 are defined in terms of the spring constant (k), pendulum lengths (Ii), masses (mi)' gravitational acceleration (g), and amplitudes (Ai) and frequency (0.) of the supports.Note that n = p = 2 in this problem.The 'differential state space fonnulation' reduces to equation ( 27) with where a ;= 02 ai' b ;= 0 2 bi, c; =. 0 2 Ci(from multiplying through by T-= 02 ), Q is the product operational matrix associated with COS(21T'T), and all submatrices are of dimension m X m.The zero imth rows of Dr, i = 1,.. .,4,correspond with the rows of Z4 and 0 that are replaced with the four initial conditions (cf.equation ( 29)), and the resulting system of 4m algebraic equations may be solved to obtain the solution.The 'differential direct fonnulation' reduces to equation ( 34) with c *l in which the four initial conditions (equation ( 36)) are inserted into the imth and (im-1)th rows of Zz and 0, i = 1,2, which correspond to the zero rows of [D If This results in a system of 2m equations which may be solved to obtain the solution.Of course, this approach becomes more efficient than the 'differential state space fonnulation' as the number m of Chebyshev tenns is increased due to the reduced number of equations to be solved.
However, this problem can be solved in a much more efficient way for large m by rewriting equations (49) as a single fourth order equation (in which pn = 4 is kept constant) and applying the 'differential direct fonnulation' which then yields only a set of m algebraic equations.From equations (49), the equivalent fourth-order equation is obtained as where QI3'Q-y,and Q&are the product operational matrices associated with cos(211''T), sin(211"T), and cos(411"T), respectively, and the equation has been multiplied through by 1611'4 -r = n4' The four initialconditionsfromequation(36)are insertedin the (m-i)throws of Z\ and 0, i = 0,...,3 corresponding to the zero rows of [DTt and the resulting system of m algebraic equations may be solved to obtain the solution.It is to be noted that, unlike the transformation to state space form which preserves the original initial conditions, a correct solution to equation (49) requires a transformation of initial conditions when using the fourth-order form (equation (52».However, although the FTM's for equations ( 49) and (52) may be different, their eigenvalues which determine the stability are the same.
For the parameter set g = 9.81 m/s2, k = 2 N/m, Ii = 1 m, mi = 1 kg, (a = 11.81 and Ci = -2.0)and Ai = 0 (no support oscillation, making the system time-invariant), the eigenvalues of the system are ::!:: 1.50*i and ::!::0.277*i and the system is neutrally stable.The motion was then analyzed for different values of the two amplitudes of support oscillation in the range Ai = -1.11m to 1.11m corresponding to bi = -10.0 to 10.0 and witha frequencyof n = 3 rad/s.For each pair of amplitudes, the FTM's for equations (49) and (52) were computed via equation (42) and their eigenvalues subsequently analyzed for stability.The results are presented in Figure 2, where 53,972 points in the (bbb2) parameter plane have been analyzed for stability.Here, the dark region represents the neutrally stable zone (when all four eigenvalues of the FTM lie on the unit circle in the complex plane) and the white region the unstable one (when at least one eigenvalue is outside the unit circle) in which destabilization occurs via the tangent route.The range on the two axes is slightly different in order to demonstrate that the oval-shaped stable region is symmetric with respect to the two 45 degree lines through the origin.This is expected since the coupled equations themselves are symmetric with identical parameter values and since, as is well known from plots similar to Figure I (cf [30] and [31]), the stability of a set of parameters is invariant to changing the sign of the periodic stiffness term.Physically this just means that the support amplitudes can always be taken as positive quantities.
In Figure 3, A2 is fixed at O.3m (b2 = 2.7) and the stability of 360,000 points in the (ajobj) parameter plane is plotted where the dark regions represent the unstable zones as in Figure 1.It is observed that the Arnold tongues which destabilize via the tangent and period doubling routes intersect the bj = 0 line (which the finite resolution does not show for some of them as in Figure 1) near aj = (q-rrITf = 9r/14,q = 1,2,...whereT = 2-rr/3 (from normalizing the {} = 3 frequency to 2-rr) which correspond to the b = 0 intersections in the single Mathieu equation.The coupling causes these points to be slightly shifted for small q, and the largest change occurs for q = 2 in which the coupling shifts the bj = 0 intersection from aj = 9.0 to aj = 10.4.However, the amount of shift quickly becomes insignificant as q increases.In addition, another route of destabiIization (due to b2 being nonzero) in which two pairs of eigenvalues of the FTM meet on the unit circle of the complex plane (but not on the real axis) and then split is exhibited in every third tongue.This destabilization route, called a Krein collision, requires at least four eigenvalues and is therefore found only in Hamiltonian systems with at least two degrees-of-freedom.Together with tangent and period doubling, these are the only three possible routes of destabilization for a linear symplectic mapping (e.g. the FfM for a linear periodic Hamiltonian system) [32].Increasing h] from 0 at a] = 11.81 in this figure corresponds to increasing h] from 0 at h2 = 2.7 in Figure 2. In either case, one travels from a stable zone to an unstable zone (due to tangent destabilization) at h] = 5.63.In Figure 4, a] is fixed at 33.0 (a2 is still 11.81, making the system asymmetric) and 21,600 points in the (hbh2) parameter plane are again analyzed for stability.Here, the dark regions are neutrally stable and the white regions unstable as in Figure 2. A ring-like region of Krein collision destabilization can be seen in the interior, while the outer region corresponds to period doubling.Increasing h] from 0 at h2 = 2.7 in this figure corresponds to increasing h] from 0 at a] = 33.0 in Figure 3.In either case, one travels from a stable zone to an unstable zone (due to Krein collision) at h] = 18.5, to another stable zone at h] = 21.5, and finally to another unstable one (due to period doubling) at h] = 39.0.
It was found that the cpu time required for the analysis of one parameter set via the fourth order 'differential direct formulation' (equation ( 53)) with a 14 Chebyshev polynomial expansion was actually more (0.37 s) than that required via the second-order 'differential direct formulation' (0.25 s).This is because the time to perform the additional matrix multiplications outweighed the time saved in solving 14 less algebraic equations.However, the cpu time required with a 50 term expansion via the fourth order form was less (0.69 s) than that required via the second-order form (1.23 s) since the time saved in solving 50 less equations outweighed the matrix multiplication time, where all times given have been averaged from 10 runs each on a SUN Sparc 20 workstation.It is expected that increasing the number n of equations would have a similar effect, and that the 'differential direct formulation' used in equivalent higher order equations (keeping pn constant) would be the fastest even for small m if the number n of equations was large.It was reported in references [15,16] that the 'integral formulation' (especially the direct approach) was several times faster than single-pass Runge-Kutta, Adams-Moulton, and Gear methods in computing the characteristic exponents of a large-order system.It is therefore to be anticipated that the 'differential formulation' would also achieve good efficiency although in the present approach pn additional linear equations are needed to be solved for an equivalent accuracy (cf section 3.2).For example, if 15 terms are used for a set of ten second-order differential equations, then in the 'differential direct formulation' one needs to solve a set of 170 algebraic equations as opposed to 150 required for the 'integral direct formulation.'The difference is undetectable in the low-order examples used in this study since a large portion of the cpu times reported above was spent in performing matrix multiplications; however, the difference would be more evident for larger systems.The initial memory overhead is exactly the same as in the 'integral formulation' since in both formulations storage of the differentiation or integration and product operational matrices is required.

CONCLUSIONS
It has been demonstrated that a set of pth order linear differential equations with periodic coefficients may be solved and analyzed for stability by expanding the periodic matrices and the solution vector in terms of shifted Chebyshev polynomials of the first kind.The use of the differentiation and product operational matrices leads to the reduction of the original set of differential equations to a set of algebraic equations in terms of the Chebyshev coefficients for the solution vector.The Hoquet Transition Matrix (FTM) may then be computed and its eigenvalues (the Hoquet multipliers) or the associated characteristic exponents analyzed for stability.Two methods were given for computing the solution vector and its derivatives.In the first, the original equations are rewritten in state space (first order) form in such a way as to completely avoid matrix inversions.The 'differential state space formulation' may then be employed which converts the problem to solving pnm linear algebraic equations in which p is the order of the equations, n is the number of equations, and m is the number of Chebyshev terms used in the expansion.The 'differential direct formulation', on the other hand, applies directly to the orginal set of differential equations and leads to solving a set of nm linear algebraic equations in terms of the unknown coefficients for the solution vector.It was shown that this approach is more efficient when either m or n is large.Unlike the similar 'integral direct formulation' used previously [15][16][17][18], it was shown that 'differential direct formulation' eliminates the need for repeated integration by parts and is much more straightforward.In both differential formulations, the special properties of the differentiation operational matrix allow the initial conditions to be included without adding additional algebraic equations.Both formulations were used in the analysis of a single Mathieu equation and of a pair of coupled Mathieu equations and stability.diagrams for certain parameter sets were constructed in which unstable zones were labeled with the type of destabilization occuring in these zones.The coupled equations were also analyzed as one fourth-order equation where the 'differential direct formulation' yielded the most efficient (for large m) method of analysis.where ar and br are Chebyshev coefficients of the functions fit) and g(t), respectively.Using equation (Ol), we can rewrite equation (02) as where Q is the product operational matrix of shifted Chebyshev polynomials of the first kind given by I am-I ::(am-Z+ am) having both a zero row and a zero column.(A matrix A is nilpotent if A k = 0 for some positive integer k.The minimum such k is called the index of nilpotence.)To find the index of nilpotence of D or DT, the following effect of the kth power of D (respectively DT) is considered.The multiplication of D (respectively DT) by itself k -1 times transforms k-1 successive subdiagonals of D (respectively superdiagonals of DT) from non-zero to zero and thus increases the number of zero rows and columns to k and decreases the rank to m -k.Since there are m -1 subdiagonals of an m X m D matrix (respectively superdiagonals of DT), the multiplication of D (respectively DT) by itself m -1 times removes all of the non-zero elements (i.e.Dm= [DTr = 0), and so the index of nilpotence is simply m.
21Twhich has period T = n and parameters a and b.Note that n = 1 and p = 2 in this problem.Transforming to state space form and normalizing with t = TT yields equation(22) with WI = 12, AI(T)that the entire equation has been multiplied by fl = 02.Followingthe procedure outlined in Section 4.1, it can be shown that the problem reduces to equation (27), where [ DT Zpn = Z2 = :;: (aI + bm

Figure 1 .
Figure 1.Neutrally stable (white) and unstable (dark) regions on the (a,b) parameter plane for the Mathieu equation where 0. = 27.The routes of destabilization consist of tangent and period doubling.