A new approach of computing Floquet transition matrix

A new approach of computing Floquet transition matrix

Computers and Structures 79 (2001) 631±635 www.elsevier.com/locate/compstruc A new approach of computing Floquet transition matrix Zhiqin Cai *, Yua...

154KB Sizes 0 Downloads 5 Views

Computers and Structures 79 (2001) 631±635

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 eciency 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 ecient, economical and viable. Ó 2001 Elsevier Science Ltd. All rights reserved. Keywords: Precise integration; Ordinary di€erential equations with periodic coecients; 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 di€erential equations with periodic coecients y_ ˆ A…t†y; …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…t†U…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 ecient 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 di€erential 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…t†y ˆ 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 ; tk‡1 †, we have f ˆ …A ÿ A0 †y ˆ

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

…5†

iˆ1

where gik , hik , rik and sik are constant vectors gik ˆ Di …yk‡1 ÿ yk †=Dt, hik ˆ Bi …yk‡1 ÿ yk †=Dt, rik ˆ ÿDi ‰yk‡1 tk =Dt ÿ yk …1 ‡ tk =Dt†Š and sik ˆ ÿBi ‰ yk‡1 tk =Dt ÿ yk …1 ‡ tk =Dt†Š. Consider a nonhomogeneous linear di€erential 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 yk‡1 ˆ T0 …Dt† yk ÿ ‰tk …aik sin itk ‡ bik cositk † ‡ cik sinitk ‡ dik cos itk Š iˆ1

m X ‡ ‰tk‡1 …aik sin itk‡1 ‡ bik cositk‡1 † ‡ cik sin itk‡1 ‡ dik cositk‡1 Š;

…7†

iˆ1

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 yk‡1 at tk‡1 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 yk‡1 is " m m X X yk‡1 ˆ In ‡ T0 wi …w1i sin itk ÿ w2i cositk †=Dt ÿ wi …w3i sin itk‡1 ÿ w4i cositk‡1 ‡ w1i sin itk‡1 =Dt iˆ1

ÿ w2i cositk‡1 =Dt†

#ÿ1 (

iˆ1 m X T0 ÿ T0 wi ‰… ÿ w1i =Dt ‡ w3i † sin itk ‡ … ÿ w4i ‡ w2i =Dt† cos itk Š iˆ1

) m X wi … ÿ w1i sin itk‡1 ‡ w2i cositk‡1 †=Dt yk ; ‡

…8†

iˆ1

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, yk‡1 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; 1ŠT . 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 ecient 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 di€erent numerical methods for the ordinary di€erential 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 di€erential 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 yˆ6 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 di€erent 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 eciency 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 ecient, economical, viable and easily applicable to any engineering problem represented by linear di€erential equations with periodic coecients.

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

References [1] Friedmann P, Hammond CE. Ecient 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 di€erential 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 di€erential equations. SIAM Rev 1984;26(3):417±21.