Seria: M E C H A N IK A z. 115 N r kol. 1230
F ranciszek A. D U L
Instytut T echniki Lotniczej i M echaniki Stosowanej Politechnika W arszaw ska
N O W A M E T O D A W Y Z N A C Z A N IA C Z Ę S T O Ś C I I P O S T A C I D R G A Ń W Ł A S N Y C H U K Ł A D Ó W F IZ Y C Z N Y C H O W IE L K IE J L IC Z B IE
ST O P N I SW O B O D Y
S treszczen ie. P rzedstaw iono now ą m eto d ę w yznaczania dużej liczby postaci w łasnych, szczególnie przydatną dla układów sym etrycznych o wielkiej liczbie stopni sw obody, N —1 0 ^ 1 0 . P okazano, iż je s t o n a znacznie efektyw niejsza od standardow ych m etod: Iteracji P odprzestrzennych i L anczosa. O bliczenia p r a ktyczne potw ierdzają przew agę now ej m etody, ro sn ącą ze w zrostem liczby stopni sw obody N. P ro p o n o w an a m eto d a m oże istotnie zm niejszyć koszt a n a l
izy m odalnej przy projektow aniu konstrukcji.
T H E N E W M E T H O D F O R C O M P U T IN G O F F R E Q U E N C IE S A N D M O D E S H A P E S O F P H Y S IC A L SY STEM S W IT H V E R Y L A R G E
N U M B E R
O F D E G R E E S O F F R E E D O M
S um m ary. T h e new conrnutational m ethod, especially su ited to th e partial solution o f veiy large, N ~ U r - f 10 generalized sym m etric eigenproblem s has b een p re s e n te d . Its effectiveness is m uch larger th a n o f th e sta n d a rd m ethods:
Subspace Iteratio n s and Lanczos. T h e results o f tests fo r various problem s confirm this advantage, w hich grows quickly with th e n u m b e r o f degrees of freed o m N. It m ay th e n b e used for considerable reducing o f cost of m odal analysis in course o f th e design process.
H O B H H M ETO II B IM H C JIEH H H U A C T O T H $0PM K.O JIEB AH H H
< M 3 IM E C K H X C H C T EM C O U E H b BOJIbUIHM 4HCJI0M C T E r iE H H C B O B O D U
r P e 3 i o M e . O n u c a H h o b h ü M e T o n p e m e H H H o u e H b 6 o n b u i H X n ~ i o 3 -h i o D . c H M M e T p t m e c K H X o 6 o 6 i u e H H H X 3 a n a u H a c o 6 c T B e H H e 3 H a U e H H 9 H B e K T O p H . ErO 3 $ $ e K T H B H 0 C T b 3 H a U H T e n b H O 6 o n b I l i a f l u e M C T a H - n a p T H H X M e T o n o B : U T e p a i t H H I l o n n p o c T p a H C T B h J I a H u o c a . I i p a K T H u e c K n e B h m H c n e H H s c o B c e M y T B e p x n a i o T s t o n p e H M y m e c T B O , K O T o p b i e p a c T e T 6 h c t p o c y B e n m t e H H e M N.
r i o s T O M y M e T o n M o x e T 6 h t b n p u n o s c e H k yM e H b iu eH H io c t o h m o c t h M o u a n b H o r o a H a n n 3 a n p H n p o e K T H p o B a H H H c j i o x h h x K O H C T p y K U H f t .
80
1. W S T Ę P
P ro jek to w an ie konstrukcji lotniczych, kosmicznych, m orskich czy budow lanych, jak rów nież inne problem y, np. zagadnienia akustyki, fizyki kw antow ej w ym agają często w yznaczenia pew nej liczby najniższych w artości własnych X; i w ektorów własnych jcj, i= l,..,M , M = 1t100, spełniających wielkie, uogólnione, sym etryczne rów nanie w łasne
K x = X M x , ( ! )
gdzie K, M są sym etrycznym i m acierzam i sztywności i bezw ładności o w ym iarze N (rów nym liczbie stopni swobody układu fizycznego), M je st nieujem nie określona.
W typow ych zadaniach m echaniki konstrukcji lub teorii pola m acierze K i M są rzadkie (liczba niezerow ych elem entów w w ierszach je st m ała, m g~ 5-r50) o raz m ają stru k tu rę pasm ow ą-^, M ijt K tj = 0 dla |i- y | > m b , nig - szerokość pasm a.
S tan d ard o w e m etody obliczeniow e, używane zwykle do rozw iązyw ania takich zadań:
M e to d a Iteracji P o dprzestrzennych (w łączona do w ielu pakietów C A D -M E S, np.
A D IN A , F E A P ) o raz M eto d a Lanczosa, są efektyw ne dla N < 102 . Je d n a k ż e d o k ład n a dyskretyzacja obiektów o skom plikow anych kształtach przy użyciu M E S-u prow adzi do zad ań o znacznie większych w ym iarach, osiągających często N -lO ^ -rlO ^ [1]. M etody stan d ard o w e są w takich przypadkach m ało efektyw ne ze w zględu na konieczność w ykonyw ania b ard zo kosztow nego rozkładu tró jk ątn eg o m acierzy
K - o M = L D L t, (2)
gdzie: L - m acierz tró jk ątn a dolna, D - m acierz diagonalna, a - liczba.
W pracy przedstaw iono nową m etodę, która nie korzysta z rozkładu macierzy i w zw ią
zk u z tym je st bardzo efektywna przy rozwiązywaniu wielkich zagadnień własnych [2].
A naliza teo rety czn a i testy wykazały jej przew agę zarów no n ad m eto d am i s ta n d a rd o wymi (w tym - n ad m eto d ą L anczosa, uw ażaną pow szechnie za najlepszą [1,3]), ja k rów nież innym i, nie korzystającym i z rozkładu macierzy. W celu zilustrow ania uniw er
salności i efektywności m etody przedstaw iono różne zagadnienia rozw iązane za jej pom ocą.
2. O P IS M E T O D Y
M e to d a po leg a na w yznaczaniu kolejnych p a r własnych p o p rz e z ¡terow anie zbioru w ektorów w dwu fazach: w stępnej, w której wyznacza się przybliżenia p a r w łasnych z m ałą dokładnością, lecz małym kosztem , i właściwej, w której kolejne w ektory w łasne
Dalsze rozważania i wnioski są słuszne również dla innych postaci macierzy.
w yznacza się z żą d a n ą , wysoką dokładnością.
Ją d re m algorytm u są iteracje odw rotne, w których kolejne przybliżenia w ektorów w łasnych w yznacza się iteracyjnie z układu rów nań
( K - a w A / ) y ("łl> = M y w , n= 0,1,... ( 3) gdzie je s t tzw. przesunięciem , a n - num erem iteracji.
M e to d a sk ład a się z dwu faz. W fazie pierwszej zbiór w ektorów (p o d p rzestrzeń ) Sg" = ] je st popraw iany iteracyjnie za p o m o cą p ro ced u r: iteracji o d w rotnych (3) z przesunięciem a ( n )*A j, R ayleigha-R itza o raz G ram a-S chm idta. Jeżeli w skaźnik stabilizacji pew nego w ektora y ^ \ zdefiniow any jak o
% y l? ) =
I
r t ? ) - \/
\ p ( y t l ) \ ,W
gdzie: n - n u m e r iteracji, Mtj*0) - iloraz R ayleigha dany w zorem
n(y) = y TK y I y TM y , (5)
spełni w aru n ek 9 { y ^ ) < r ( r - p róg przełączania), to przechodzi się do fazy drugiej, w której za p o m o cą iteracji odw rotnych (3) z przesunięciam i a^n ^=#x(yl"*) oblicza się przybliżenie w ek to ra w łasnego z ż ąd an ą dokładnością c
II e(y(? ) || = || K y (? - n ( y ? ) M y\? || < e . (6) O bliczony w e k to r własny zastępow any je st w zbiorze S g 1 nowym w ek to rem s ta rto wym i cały p ro ces je s t pow tarzany aż wyznaczone zo stan ą wszystkie poszukiw ane pary w łasne.
Iteracje o d w ro tn e w ym agają rozw iązania układu (3) z m acierzą n ieokreśloną i praw ie osobliw ą d la -*■ Aj. Aby uniknąć rozkładu tró jk ątn eg o (2) do jego ro z
w iązania używa się specjalnych p ro ced u r iteracyjnych [4]. Isto tn ą p o p raw ę efektywności uzyskuje się sto su jąc tzw. skalow anie (preconditioning). Szczegóły przedstaw ione są w artykule [2],
M e to d a zrealizow ana je st w postaci autonom icznego p ak ietu , d o którego należy dołączyć trzy procedury, obliczające w ektory a := K x, b:= M y, c:= P "3z (P jest m acierzą skalującą), dla w ektorów x, y i z przekazyw anych przez p ro g ram wywołujący. N a w ej
ściu należy po d ać: liczbę poszukiwanych postaci własnych M i do k ład n o ść c.
M e to d a je s t b a rd z o dobrze dostosow ana do obliczeń rów noległych, gdyż w obu jej fazach wykonywanych je st szereg niezależnych procesów obliczeniow ych.
82
3. Z A S T O S O W A N IE P R A K T Y C Z N E M E T O D Y
M e to d a um ożliwia rozw iązyw anie dowolnych, sym etrycznych zagadnień własnych.
Je d n a k jej przew aga n ad innymi m etodam i uw idacznia się dla zadań, w których N > 1 0 0 0 i m ^ N ^ - r N ® ' ^ . W arunki te są zazwyczaj spełnione w zagadnieniach teorii pola (np. akustycznego) o raz m echaniki konstrukcji, zw łaszcza trójw ym iarowych.
M e to d a um ożliw ia w yznaczenie sztywnych stopni swobody (np. sam olotów lub w iel
kich konstrukcji orbitalnych), ja k rów nież postaci drgań maszyn w irnikowych (wirniki śmigłowców, siłowni w iatrowych, turbin gazowych, silników odrzutow ych, etc.). M odel obiek tu m oże uw zględniać tzw. tłum ienie proporcjonalne,
F D( q ) - C q(t) , C = ctK + PAf , a , p - współczynniki. (7) R ozw iązyw ane m ogą być też zadania z osobliwą m acierzą M , pow stającą przy pom ijaniu w m odelu fizycznym niektórych efektów bezw ładnościow ych. Innym i obszaram i zastosow ań m etody są: zadania pełne, tzn. w yznaczenie wszystkich M = N w ektorów w łasnych o raz zagadnienia słabo niesym etryczne [2]. Z a p o m o cą opisanej m etody rozw iązano w iele zadań, dotyczących głów nie układów m echanicznych;
w yznaczono m.in.:
- P o stacie akustyczne w nętrz trójwymiarowych kabin sam ochodów , m odelow anych liniowymi i kw adratow ym i elem entam i skończonymi; N =956-r20186; M=204-100;
- Postacie drgań płyt m odelow anych elem entam i herm itow skim i; N =100-r6000, M —80.
- P ostacie d rg ań trójw ymiarowych konstrukcji belkowych; N = 720-1-10368, M —50.
- Postacie drgań służące do analizy dynamicznej konstrukcji złożonych z w ielkich u kładów liniowych połączonych małymi układam i o nieliniowych charakterystykach sprężystych i tłum iących; N = 1120, 2460, 3847; M =10-^35.
Rozw iązyw ano także inne problem y, na przykład:
- Z ag ad n ien ia chem ii kw antow ej: wielkoskalow e oddziaływ ania funkcji falowych atom ów i m olekuł; N = 10004-40000, M =10-r350.
- P ełn e zag ad n ien ia w łasne, np. N = M = 7 2 9 ;
- E k strem aln ie wielkie zagadnienia w łasne: N =103823, M = 3 ; N =216000, M = 2 .
4. P O R Ó W N A N IE Z IN N Y M I M E T O D A M I
P rzydatność praktyczna każdej m etody zależy zarów no od jej efektyw ności, ja k też w łasności użytkowych: uniw ersalności oraz niezawodności. M etody w ykorzystujące ro z
kład tró jk ątn y (2) (zw łaszcza m e to d a L anczosa) są b ard zo efektywne, gdy, N < 1 0 0 0 lub
szerokość p asm a m b je st m ała. N ie są je d n a k efektyw ne dla dużych N i m b, gdyż ich koszt (m ierzony liczbą operacji), rzędu ~ N 2-2-rN3, rośnie dużo szybciej niż w przedsta
wionej m etodzie, gdzie je st on rzędu - N ^ - r N 1'6 . P o n ad to m acierz L zajm uje N m b jed n o ste k pam ięci, podczas gdy K i M zajm ują ich łącznie 2NmQ. D la typowych w artości: m g = 1 0 i m b = 200 m acierz L zajm uje 10 razy więcej pam ięci niż K i M.
S tosunek ten rośnie ze w zrostem w ymiaru N, gdyż m b ~N® ‘3-rN® , podczas gdy nig—
const. M eto d a L anczosa w ymaga dodatkow o pam ięci ~50N -r300N dla w ektorów ro b o czych.
Przew agę p rzedstaw ionej m etody dla zadań o dużej liczbie stopni swobody całko
wicie potw ierdzają obliczenia praktyczne (R ys.l): ju ż dla N = 4000, M = 1 0 jest ona około dziesięć razy szybsza niż Iteracje P odprzestrzenne i trzy razy szybsza niż m etoda Lanczosa. Pierw sze postacie w yznaczane są w czasie stanow iącym u łam ek czasu niez
b ęd n eg o do ro zk ład u m acierzy (np. 1/100 dla N = 10000).
N iezaw odność m etody, po d o b n ie ja k Iteracji P odprzestrzennych, je s t wysoka, w przeciw ieństw ie do m etody L anczosa, któ ra nastręcza kłopotów związanych z w ielokrotnym i w artościam i własnymi, często spotykanym i przy analizie konstrukcji sym etrycznych [2,3],
W yniki testów pokazały, że inne m etody, nie wykorzystujące ro zk ład u (2), np.
T R A C M N , S IR Q IT , S R Q M C G , nie są tak efektywne, p oza tym są często wrażliwe na rodzaj zad an ia (m e to d a D avidsona) lub przybliżenia początkow e (T L IM E ) [2],
R y s.l . P o ró w n an ie kosztów względnych w yznaczania postaci w łasnych trójw ym iaro
wych konstrukcji belkowych. P okazano też koszt rozkładu m acierzy
F ig .l . C o m p ariso n o f relative costs o f finding o f m o d e sh ap es fo r 3-D building- type p ro b lem s. T h e cost o f factorization is also shown
84
5. P O D S U M O W A N IE
Z e w zględu n a fakt, iż projektow anie wielu obiektów technicznych w ymaga rozw iązyw ania coraz większych zagadnień własnych, m ożna spodziew ać się, że p rzed staw io n a m e to d a znajdzie zastosow anie w system ach C A D , zważywszy jej zalety:
dużo niższy (często 10-krotny) koszt obliczeń w porów naniu z m etodam i standardow ym i, uniw ersalność, niezaw odność o raz łatw ość użycia.
L IT E R A T U R A
[1] M atth ies H .G .: A Subspace Lanczos M eth o d for th e G en eralized Sym m etric E ig en p ro b lem , C o m p u ters & Structures, vol.21, no.1/2, 1984, ss.319-325.
[2] D ul F.A ., A rczew ski K.: T h e T w o-Phase M ethod for Finding a G re a t N u m b e r of E ig en p airs o f th e Sym m etric o r W eakly N on-Sym m etric L arg e E igenvalue p ro b lem s, Jo u rn al o f C o m p u tatio n al Physics (w druku, 1993).
[3] P a rle tt B.N.: T h e Softw are Scene in th e E xtraction o f E igenvalues from S parse M atrices, SIA M Jo u rn al on Scientific an d Statistical C om p u tatio n s, vol.5, n o .3, 1984, ss.590-604.
[4] P aig e C.C., S aunders M .A.: Solution o f Sparse Indefinite System s o f L in ear E q u atio n s, SIA M Jo u rn al on N um erical Analysis, v ol.12, no.4, 1975, s.617-629.
R ecen zen t: D r hab. inż. Jerzy Świder W płynęło do R edakcji w grudniu 1993 r.
A b stract
T h e designing o f space structures, aircraft, ships, buildings, oil pipes, etc., as well as o th e r p ro b lem s o f physics, re q u ire often com puting a few eigenpairs Xj, x-v i= l,..,M , M = 1-rlOO, o f large (N —lO^-rlO ), generalized sym m etric e ig en p ro b lem (1). T he sta n d a rd co m p u tatio n al m ethods: Subspace Iteratio n s an d Lanczos, perfo rm facto rizatio n (2), th u s a re n o t efficient fo r large p roblem s b ecau se o f huge m em ory re-
i r i /'/'
1
q u irem en ts ( ~ N -rN ) and unacceptable cost o f co m putations ( ~ N -rN ).
The new m eth o d presented here, being factorization-free, is especially suited fo r solving very large eigenproblems, du e to m oderate cost ( —N ^ - r N ), sm all m em ory requirements ( ~ N ) a n d high level o f reliability [2],
T h e m e th o d consists o f tw o phases: in th e first phase th e set Sjj.'0 = ]» >s ite ra te d using: Inverse Iteratio n s (3) w ith co n stan t shift
< /n Q X j, R ayleigh Q u o tien t and G ram -Schm idt. A fter k-th vector is stabilized, ( K y ^ )< t
(4), th e second p h ase is p erfo rm ed using Inverse Iteratio n s with <7 ■n )=/r(_y(J")) (5) to o b tain the final accuracy lk(y^)ll< e (6). T he process is th en re p e a te d with new starting vector until all M req u ired eigenpairs are com puted. T h e indefinite, alm ost singular system s (3) are solved iteratively [4] with preconditioning [2], T he m eth o d is well suited to th e p arallel com p u tatio n s.
T h e m eth o d is g en eral p urpose, b u t it is especially suited to p ro b lem s with N > 1 0 0 0 and large b an d w id th (- N ^ -5 - N * * '^ ). It is able to co m p u te rigid m odes o f fre e con
structions as well as th ese o f rotating m achinery. T h e p ro p o rtio n al dam ping (7) or singular m ass m atrices m ay b e taken into account. Som e weakly non-sym m etric eigen- p ro b lem s can be solved [2], H e re is the collection o f problem s solved by the pro p o sed m ethod:
- A coustics m o d es o f th e 3-D interiors discretized by FE , N = 956-5-20186, M = 20-5-100;
- M o d e sh ap es o f p lates discretized by herm itian F E ; N = 100-5-6000, M —80.
- V ib ratio n m o d es o f 3-D buildings com posed o f beam s: N=720-5-10368, M —50.
- M o d e sh ap es fo r dynam ic analysis o f mixed stru ctu res (larg e-lin ear jo in t by sm all n on -lin ear); N = 1 1 2 0 , 2460, 3847; M =1 0-5-3 5;
- P ro b lem s o f q u an tu m chemistry: the large scale configuration interactions of electro n ic w avefunctions o f atom s and m olecules; N = 1000-5-40000, M=10-5-350;
- Full (e.g., N = M = 7 2 9 ) and extrem ely large (e.g., N = 216000, M = 2 ) eigenproblem s.
T h e effectiveness o f the m ethod is very high, its superiority grow s quickly w ith the dim ension N (F ig .l), e.g., for N =4000, M = 1 0 it is a b o u t ten tim es faster th an the Subspace Ite ra tio n s and th re e tim es faster th a n Lanczos. T h e first eigenm ode is co m p u ted w ith th e cost being a sm all fraction o f th a t o f factorization. T h erefo re, it can be u sed fo r significant reduction of the cost o f m odal analysis, an im p o rtan t p a rt o f the design process.
