- Email: [email protected]

www.elsevier.com/locate/compstruc

A new approach of computing Floquet transition matrix Zhiqin Cai *, Yuanxian Gu, Wanxie Zhong Research Institute of Engineering Mechanics, Dalian University of Technology, Dalian 116 024, People's Republic of China Received 2 December 1999; accepted 11 May 2000

Abstract A numerical solution based on the precise integration method is presented for computing the Floquet transition matrix (FTM), and dealing with the stability of the periodic system. The numerical properties of this solution are illustrated by comparing the numerical results and eciency of two methods: (1) Hamming's predictor±corrector, (2) Runge±Kutta. It is concluded that the precise integration method in single-pass is the most ecient, economical and viable. Ó 2001 Elsevier Science Ltd. All rights reserved. Keywords: Precise integration; Ordinary dierential equations with periodic coecients; Floguet transition matrix; Stability problem

1. Introduction The equations of motion of a wide class of dynamical systems including rotorcraft [1] are formulated as the ordinary dierential equations with periodic coecients y_ A ty; 1 where y is the unknown n-dimensional vector, A the n n matrix, A t A t T , y_ dy=dt. Based on the Floquet±Liapunov theorem [2], the stability of such a system is analyzed by ®nding the most positive real part of the eigenvalues of the Floquet transition matrix (FTM). FTM is a n n transition matrix U t at the end of one period, and satis®es _ U t A tU t; U 0 In : 2 For the FTM, there were principally two classes of computational approaches. In the ®rst class of approaches, Hsu [3] has developed various methods for approximating the transition matrix during one period; the most ecient one consists of approximating the periodic matrix A(t) by series of step functions. In the second class of approaches, the direct numerical integration for the ordinary dierential equations is directly applied; the Hamming's fourth order predictor±corrector method is the most economical and is more economical than Hsu's methods [4]. In this paper, precise integration method [5] is extended for computing FTM and compared with the above two classes of computational approaches. 2. Application of the precise integration method It is convenient to rewrite Eq. (1) in the form y_ A ty A0 y f ; f A ÿ A0 y; P where A ÿ A0 Di sin it Bi cos it, and A0 , Di , Bi i 1; 2; . . . ; m are the constant matrices. *

Corresponding author.

0045-7949/01/$ - see front matter Ó 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 0 0 ) 0 0 1 6 9 - 3

3

632

Z. Cai et al. / Computers and Structures 79 (2001) 631±635

The step-by-step integration is necessary. Let the length of a time step be Dt T =M, the solution vectors yk at time t0 0, t1 Dt; . . . ; tk kDt; . . . ; tM T are to be calculated. First, matrix T0 is calculated with the precise integration method [5] T0 W Dt exp A0 Dt;

4

_ where W t A0 W t, W 0 In (identity matrix), and T0 can be considered to be precise. Next, approximating that y in f of Eq. (3) is linear within a time step tk ; tk1 , we have f A ÿ A0 y

m X gik t sin it hik t cosit rik sin it sik cos it;

5

i1

where gik , hik , rik and sik are constant vectors gik Di yk1 ÿ yk =Dt, hik Bi yk1 ÿ yk =Dt, rik ÿDi yk1 tk =Dt ÿ yk 1 tk =Dt and sik ÿBi yk1 tk =Dt ÿ yk 1 tk =Dt. Consider a nonhomogeneous linear dierential equation y_ A0 y t r1 sin it r2 cosit r3 sin it r4 cosit;

6

where an integer i 6 0 and rj j 1; . . . ; 4 are constant vectors. Its special integral is in the form t a sin it b cos it c sin it d cosit, a wi ÿA0 r1 =i r2 , b ÿwi A0 r2 =i r1 , c wi A0 a=i ÿ A0 r3 =i r4 ÿ b, and d ÿwi A0 r4 =i ÿ A0 b=i r3 ÿ a, wi iIn A20 =iÿ1 . Then, integrating Eqs. (3) and (5) gives ( ) m X yk1 T0 Dt yk ÿ tk aik sin itk bik cositk cik sinitk dik cos itk i1

m X tk1 aik sin itk1 bik cositk1 cik sin itk1 dik cositk1 ;

7

i1

where aik wi ÿA0 gik =i hik , bik ÿwi A0 hik =i gik , cik wi A0 aik =i ÿ A0 rik =i sik ÿ bik , and dik ÿwi A0 sik = i ÿ A0 bik =i rik ÿ aik . When the vector yk1 at tk1 is calculated from Eqs. (5) and (7), gik , hik , rik and sik in Eq. (5) should be substituted into Eq. (7), and the resulting equation involving only yk1 is " m m X X yk1 In T0 wi w1i sin itk ÿ w2i cositk =Dt ÿ wi w3i sin itk1 ÿ w4i cositk1 w1i sin itk1 =Dt i1

ÿ w2i cositk1 =Dt

#ÿ1 (

i1 m X T0 ÿ T0 wi ÿ w1i =Dt w3i sin itk ÿ w4i w2i =Dt cos itk i1

) m X wi ÿ w1i sin itk1 w2i cositk1 =Dt yk ;

8

i1

where w1i wi A0 w3i =i w4i , w2i wi A0 w4i =i ÿ w3i ; w3i ÿA0 Di =i Bi ; Eq. (5) gives gik , hik , rik and sik . Finally, yk1 can be obtained from Eq. (7).

w4i A0 Bi =i Di . Substituting yk 1 into

3. N-pass and single-pass schemes FTM is the n n matrix solution of Eqs. (1) or (2), and thus there are two computational approaches. The ®rst approach is referred here as the N-pass scheme. In this approach, Eq. (1) is solved n times for n initial states which comprise the n columns of identity matrix. The second approach is called the single-pass scheme. In this approach, the n n FTM as a n2 1 vector is computed only once by integrating Eq. (2) with the initial state 1; 0; 0; . . . ; 0; 0; 1; 0; 0; . . . ; 0; . . . ; 0; . . . ; 0; 0; 1T . The machine time saving through the single-pass scheme is shown [4]. It is seen that the fourth order Runge±Kutta method in single-pass is more ecient than Hsu's step function method [1]. It is seen that the detail comparisons of CPU time by N-pass and single-pass schemes, using the dierent numerical methods for the ordinary dierential equation system, are shown [4]. They demonstrate that HammingÕs fourth order predictor±corrector method in single-pass is most economical with respect to three signi®cant ®gure accuracies. Accordingly, the numerical properties of the precise integration method are illustrated by comparing the

Z. Cai et al. / Computers and Structures 79 (2001) 631±635

633

computational results and CPU time by N-pass and single-pass schemes with HammingÕs fourth order predictor± corrector method and the fourth order Runge±Kutta method up to the same accuracy. 4. Numerical examples As two examples, consider two dierential equations with analytic solutions [6] First example: ÿ2t e cost ÿ sin t ÿ2 cos2 t ÿ1 ÿ sin 2t y y; T p; its solution U t ; 1 ÿ sin 2t ÿ2 sin2 t eÿ2t sin t cost

U 0 I2 :

Fig. 1. (a)±(f) y(i), ith component of the solution vectors, ANALYTIC, analytic solution, HAMMING, solution Hamming method, RUNGE±KUTTA, solution with Runge±Kutta method, PIM, solution with the precies integration method.

634

Z. Cai et al. / Computers and Structures 79 (2001) 631±635

Second example: 2 ÿ2 cos2 t 6 1 ÿ sin 2t y6 4 3 0

ÿ1 ÿ sin 2t ÿ2 sin2 t 0 3

where the matrix 2

eÿ2t cost 6 eÿ2t sin t U t 6 4 u 3eÿt cost 1 u2 3eÿt sin t

0 0 ÿ3 2 cos2 t 1 sin 2t

ÿ sin t cost 1 eÿ3t sin t 1 ÿ eÿ3t cost

3 0 7 0 7; ÿ1 sin 2t 5 2 ÿ3 2 sin t

0 0 eÿt cos t eÿt sin t

0 0

T p

3

7 7; ÿeÿ3t sin t 5 ÿeÿ3t cost

U 0 I4

is a fundamental solution, and u1 ÿ3eÿ2t cos t, u2 ÿ3eÿ2t sin t. In order to compare dierent methods with respect to the same accuracy and to achieve the comparable program execution time, all the algorithms have built-in mechanisms of altering the step-size in response to the stipulated tolerance. For the given step-size Dt and tolerance e, when the error exceeds this given value, the program halves the interval until an acceptable error level is reached. The computational results are shown in Fig. 1(a), (b), (d) and (e), where the certain components of solution vectors, at step-size Dt T =40 0:0785 and T =20 0:1571, are plotted. The computed values are shown in Fig. 1(c) and (f), using the automatic step-size control up to the stipulated tolerance e 1:0 10ÿ5 and corresponding to the known conditions in Fig. 1(b) and (e), respectively. From Fig. 1(a), (b), (d) and (e), it is seen that the computational errors of Hamming's predictor±corrector method are smaller than ones of Runge±Kutta method, which is identical with the conclusion [4]. The errors of the precise integration method are the smallest. The three numerical methods with built-in step-size control are taken up again in Tables 1 and 2, respectively, which include the execution time for computing FTM by N-pass and single-pass with respect to the same tolerance e 1:0 10ÿ5 , at T/40, T/20, T/10, T/5 and T/3. As seen from Tables 1 and 2, for both N-pass and single-pass, CPU time with the precise integration method takes the least amount, and Runge±Kutta method, the highest. When the selected step-size is larger (for example Dt T =3) the saving in the precise integration method through single-pass is the most and close to 1=36 in the Hamming's method, whereas there is no computational result using Runge±Kutta method at this step-size. Further comparisons demonstrate two points. (1) As expected, for the precise integration method in single-pass, the saving in machine time Table 1 CPU time seconds on Pentium 100 for computing FTM of example 1 Step size

N-pass

Single

e1

e2

Total time

CPU time

PIM

0.07854 0.15708 0.31416 0.62832 1.0472

0.16 0.16 0.11 0.16 0.27

0.16 0.16 0.11 0.22 0.28

0.32 0.32 0.22 0.38 0.55

0.17 0.16 0.16 0.22 0.33

HAMM

0.07854 0.15708 0.31416 0.62832 1.0472

0.22 0.27 0.55 0.99 1.76

0.55 0.88 1.64 3.13 5.22

0.77 1.15 2.19 4.12 6.98

1.15 1.98 3.74 6.92 12.08

RKM

0.07854 0.15708 0.31416 0.62832 1.0472

0.44 0.77 1.32 ± ±

1.54 3.13 6.2 ± ±

1.98 3.9 7.52 ± ±

3.4 6.64 13.23 ± ±

Note: ± abnormal program termination; RKM, Runge±Kutta method; PIM, the precise integration method; HAMM, Hamming method.

Z. Cai et al. / Computers and Structures 79 (2001) 631±635

635

Table 2 CPU time seconds on Pentium 100 for computing FTM of example 2 Step size

N-pass e1

Single e3

e4

PIM

0.07854 0.15708 0.31416 0.62832 1.0472

0.93 0.66 0.6 1.31 2.42

0.93 0.66 0.61 0.88 1.38

0.88 0.55 0.43 0.55 0.76

0.93 0.66 0.6 0.77 1.31

3.67 2.53 2.24 3.51 5.87

0.94 0.66 0.66 1.43 2.58

HAMM

0.07854 0.15708 0.31416 0.62832 1.0472

1.15 1.98 3.57 6.81 11.53

1.48 2.53 4.67 8.9 13.78

0.72 1.05 1.81 3.24 5.11

0.5 0.83 1.64 3.29 5.83

3.85 6.39 11.69 22.24 36.25

6.1 10.82 20.21 38.84 63

RKM

0.07854 0.15708 0.31416 0.62832 1.0472

3.13 5.66 8.07

4.29 8.07 14.39 ± ±

1.59 2.75 4.56 5.22 ±

1.42 2.58 3.9 ± ±

10.43 19.06 30.92 ± ±

18.35 34.72 62.12 ± ±

± ±

e2

Total time

CPU time

Note: ± abnormal program termination; RKM, Runge±Kutta method; PIM, the precise integration method; HAMM, Hamming method.

substantially increases with increasing n and is equal to about n-fold in N-pass. (2) Rather expected, for the Hamming and Runge±Kutta method, the execution time by single-pass is more than the total time by N-pass with respect to e 1:0 10ÿ5 , and the larger the step-size, the more is the excess of CPU time by single-pass over N-pass.

5. Conclusion The numerical results are illustrated by comparing the eciency of two methods: (1) Hamming's predictor±corrector, (2) Runge±Kutta. It is concluded that the CPU time for computing FTM with the precise integration method in single-pass is the least and is about 1=n in N-pass with respect to the same accuracy. That demonstrates that the precise integration algorithm in single-pass is the most ecient, economical, viable and easily applicable to any engineering problem represented by linear dierential equations with periodic coecients.

Acknowledgements This work was supported by the National Natural Science Foundation of China (19 990 510).

References [1] Friedmann P, Hammond CE. Ecient numerical treatment of periodic systems with application to stability problems. Int J Numer Meth Engng 1977;11:1117±76. [2] Coddington EA, et al. Theory of ordinary dierential equations. New York: McGraw-Hill, 1955. [3] Hsu CS. On approximating a general linear periodic system. J Math Anal Appl 1974;45:234±51. [4] Gaonkar GH, Simha Prasad DS, Sastry D. On computing Floquet transition matrices of rotorcraft. J Am Helicopter Soc 1981;26(3):56±61. [5] Zhong WX, Zhu J, Zhong XX. A precise time integration algorithm for nonlinear systems. Proc, Third World Congr Computat Mech, Japan. 1994. Vol. 1 (pp. 12±17). [6] Arthur de Kleine H. A note on the asymptotic stability of periodic solutions of autonomous dierential equations. SIAM Rev 1984;26(3):417±21.