CURS 4 - Partea II
-
Upload
zuzu-c-bubu -
Category
Documents
-
view
227 -
download
3
description
Transcript of CURS 4 - Partea II
-
a.ch. Noiembrie 2008
1
CURS 4 METODE NUMERICE PENTRU PROBLEMA DE VALORI PROPRII -----------------------------------------------------------------------------------------------------------
Partea II 1. Algoritmul QR
2. Problema generalizat de valori proprii.
1 Algoritmul QR 1.1 Metoda
QR este unul din cei mai utilizai algoritmi pentru a determina toate valorile proprii ale unei matrici. n numele metodei, Q desemneaz o matrice ortogonal, iar R o matrice superior-triunghiular (R vine de la right).
Considerm o matrice real nesingular A, simetric sau nu.
Esena algoritmului QR const n urmtoarea proprietate:
Dac A este factorizat n produsul QRA = , unde Q este nesingular, atunci matricea produsului n ordine invers, RQA = , are aceleai valori proprii ca i A (fiind similar cu A)
ntr-adevr, avem: AQR 1= , i AQQA 1= .
Algoritmul QR realizeaz factorizri succesive ale irului de matrici { kA }, unde AA =1 , definite de:
111 RQA = ; 112 QRA = ;
222 RQA = ; 223 QRA = ;
kkk RQA = ; kkk QRA =+1 ;
-
a.ch. Noiembrie 2008
2
Matricile kA , i kk RQ , , au urmtoarele proprieti:
1) Toate kA au aceleai valori proprii ca i A.
2) Dac A este simetric, sau tridiagonal, sau A are forma Hessenberg superioar, matricile kA vor pstra aceste forme.
3) Fie kk QQQQ K21= , atunci (prin inducie): kkk QAQA 111 + = .
4) Fie 1RRR Kkk = , atunci: kkk )( 1ARQ =
Matrici Hessenberg i Househooder
Forma Hessenberg (superioar) este o matrice avnd elementele zero sub sub-diagonal. De exemplu, o matrice Hessenberg 66 arat astfel:
=
6665
5554
4443
3332
2221
1611
0000000
000
aa
aa
aa
aa
aa
aa
K
K
MK
K
A
O matrice Householder nn este definit prin
TwwIP 2= ,
unde I este matricea unitate nn, iar w este un vector cu n coordonate, normalizat
relativ la norma euclidian, adic 1|||| 2 =w , sau 1=wwT , i altfel arbitrar.
Explicit,
-
a.ch. Noiembrie 2008
3
][210
01
1
1
nj
n
i www
w
w
w
KK
M
M
K
MOM
K
=P
jiijij wwp 2=
Matricea P este simetric i este propria sa invers (este unitar), ntruct avem
PP =T , i PP =1 .
Dac se consider un vector arbitrar v, definim
2||||/ vvw = , i P devine:
22||||
2v
vvIPT
= (1)
Sau, notnd
>
-
a.ch. Noiembrie 2008
4
)1(2|||| exxv m= .
Rezult
)1(2|||| exPx m=
(Proprietatea se demonstreaz astfel: se calcueaz >< xvxv , i se arat c avem >< xvvv ,2, ; se calculeaz apoi, Px.)
1.2 Algoritm:
Algoritmul parcurge urmtoarele etape:
I. Reducerea lui A la forma tridiagonal dac A este simetric, sau la forma Hessenberg dac A este nesimetric.
II. Reducerea lui A la forma triunghiular, sau bloc- triunghiular (v. mai jos). III. Descompunerea QR.
Algoritmul QR propriu-zis, const n Etapele II i III.
a) Etapa I se realizeaz nainte de algoritmul QR pentru a aduce matricea A la o forma care reduce efortul de calcul. Aceasta, pentru c numrul de operaii pe un pas k al iteraiei QR (Etapele II i III), aplicat unei matrici nn este de ordinul (Ralston & Rabinowitz, 1983):
3n pentru o matrice general;
2n pentru o matrice Hessenberg;
n pentru o matrice tridiagonal.
Etapa I se realizeaz prin pre- i post-multiplicare cu matrici Householder kP .
b) Etapa II se realizeaz prin premultiplicare cu matrici Householder kP , cum urmeaz:
APA 11 = ; APPAPA 12122 == ; ; APPPA 1211 K = nn .
-
a.ch. Noiembrie 2008
5
Matricea 1nA are forma triunghiular i este matricea R.
ntr-adevr: S notm 121 == nPPPQ K , avem AQA Tn =1 . Q este o matrice ortogonal, deci avem AQA =
1n . Cum Q este ortogonal i 1nA este superior triunghiular, urmeaz c 1= nAR .
Atunci, Etapele II i III se realizeaz prin urmtorii pai de iterare:
1) Pune AA =)1(
- Triangularizeaz )1(A :
APA )1(1)1(1 = ; )1(1)1(2)1(2 APA = ; ; )1( 2)1( 1)1( 1 = nnn APA
- La sfritul pasului, avem:
)1(1
)1(
= nAR ; )1(1
)1(2
)1(1
)1(
= nPPPQ K .
2) Calculeaz )1()1()2( QRA =
- Triangularizeaz )2(A :
)2()2(1
)2(1 APA = ; )2(1)2(2)2(2 APA = ; ; )2( 2)2( 1)2( 1 = nnn APA
- Pune:
)2(1
)2(
= nAR ; )2(1
)2(2
)2(1
)2(
= nPPPQ K .
...
k) Calculeaz )1()1()( = kkk QRA )2( k
- Triangularizeaz )(kA :
)()(1
)(1
kkk APA = ; )(1)(2)(2 kkk APA = ; ; )( 2)( 1)( 1 knknkn = APA
- Pune:
)(1
)( kn
k
= AR ; )( 1)(
2)(
1)( k
n
kkk
= PPPQ K
-
a.ch. Noiembrie 2008
6
...
Calculaia continu pn cnd elementele triunghiului sub-diagonal n )( 1knA sunt
aproximativ zero. Atunci, valorile proprii sunt pe diagonala lui )( 1knA .
Mai precis, avem:
1) Dac A are valori proprii reale, atunci: A se reduce la forma triunghiular. (calculaia se termin cnd toate elementele sub-diagonale sunt aproximativ zero).
n particular, dac A este simetric, matricea se reduce la forma diagonal. (toate elementele non-diagonale vor fi aproximativ zero). Reamintim c o matrice simetric are valori proprii reale.
2) Dac A (nesimetric) are o pereche de valori proprii complexe (conjugate), pe diagonal rmn submatrici blocuri 22, ale cror valori proprii sunt perechea de valori proprii conjugate ale lui A. 3) Dac A (nesimetric) are valori proprii i i || i este de multiplicitate m,
atunci )(kA are un bloc diagonal de ordinul m, ale crui valori proprii sunt i .
(V. Wilkinson (1965), Ralston & Rabinowitz (1983), Meglicki (2001).)
Deci, cu excepia cazului 3 cu 2>m ( 2>m valori proprii distincte i de module egale), pentru k suficient de mare, matricea )(kA are urmtoarea form:
=
s
r
k
B
B
C
A
KK
OMMMM
OMM
000
0
1
2
1
)(
,
-
a.ch. Noiembrie 2008
7
unde
= )(
22)(
21
)(12
)(11
jj
jj
j bbbbB , iar C conine elementele superioare. Avem nsr =+ 2 .
Valorile i i blocurile jB , pot apare n orice ordine pe diagonal.
Not
Cazul 3 se poate evita n practic, prin algoritmul QR cu deplasare v. mai jos
1.3 Algoritmul QR cu deplasare:
Presupunem c A are valori proprii de module distincte 0|||||| 21 >>>> n K . Se art c viteza de convergen la zero a elementelor sub-diagonale, i cea de convergen la valorile proprii a elementelor diagonale, depind de rapoartele
kii )/( 1 + , 11 ni .
Dac dou valori proprii i i 1+i sunt apropiate una de alta, convergena va fi nceat. Atunci, se utilizeaz urmtoarea tehnic pentru a accelera convergena. Se
aplic o deplasare ks la valorile proprii (acestea devin ki s ), la fiecare etap k. Exemplu: 1.2,2.2 1 == +ii ; avem 95.0/1 =+ ii . Cu 0.2=s , raportul are
valoarea 5.02.0/1.0 = .
Adic, punem
AA =1 , i
)()()( kkk
k s RQIA = ,
IQRA kkkk s+=+ )()()1(
Exist dou strategii pentru a alege deplasarea ks ( ks est o aproximaie a lui n ):
1) )(knnk as =
Pentru o matrice simetric (real), aceast strategie asigur o convergen cubic, chiar n prezena unor valori proprii multiple (Wilkinson, 1965).
2) Se calculeaz valorile proprii ale submatricii 22 cea mai de jos din )(kA :
-
a.ch. Noiembrie 2008
8
)()(1,
)(,1
)(1,1
knn
knn
knn
knn
aa
aa
Dac valorile proprii sunt reale, se alege:
ks = valoarea proprie care este cea mai apropiat de )(knna .
Dac valorile proprii sunt complexe, se ia:
ks = partea real a valorii proprii.
3) Experimentele numerice au condus la o a treia strategie, i anume:
nias kiik ,1),min( )( ==
Pentru unele matrici simetrice, aceast strategie se dovedete cea mai bun,
conducnd la un numr mai mic de iteraii, i o eroare mai mic n Ay y.
1.4 Programul QR Programul QR din ANA, calculeaz sistemul propriu (valori i vectori proprii), n succesiunea urmtoare:
I. Transformarea preliminar a matricii A (Opional).
II. Valorile proprii, prin reducerea matricii la forma bloc-triunghiular. (Programul recunoate numai blocuri 22 .)
III. Vectorii proprii (Opional).
Se calculeaz, mai nti, vectorii proprii ai matricii bloc-triunghiulare. Acetia se aduc, prin transformarea de similaritate, la vectorii proprii ai lui A.
IV. Rafinarea sistemului propriu prin iteraie invers (Opional).
V. Verificarea sistemului propriu (Opional):
Precizia sistemului propriu se verific prin listarea erorii maxime (n modul), n yAy (unde este valoarea proprie, iar y vectorul propriu asociat cu ).
-
a.ch. Noiembrie 2008
9
Exemplu 1: Matrice simetric Considerm matricea:
=
12612332
631031231
A
Programul QR, rulat cu eps = 1E-7 (toleran pentru elementele non-diagonale), eps_lambda = 1E-6 (tolerana la ), deplasare strategia 1, produce urmtoarele rezultate (matricile sunt listate n format 1pG15.7, iar valorile proprii n format 1pG24.16):
Matrix A - Tridiagonal form:
1.000000 3.741657 -4.0811097E-16 2.8277989E-16
3.741657 2.785714 -5.246476 0.000000
-4.0811097E-16 -5.246476 10.19927 -4.479586
2.8277989E-16 0.000000 -4.479586 1.015014
Iteration 19:
Reduced Matrix - Block-Triangular Form
14.32951 8.3206414E-08 -1.0862088E-10 6.0791490E-16
8.3206414E-08 4.456957 -3.9733936E-03 -2.9646362E-11
-1.0862122E-10 -3.9733936E-03 -3.415088 -2.2709916E-08
-8.1045051E-19 -2.9646498E-11 -2.2709916E-08 -0.3713752
Eigenvalues
1 14.32950642539378
2 -3.415090280621962
3 4.456959098788067
4 -0.3713752435599103
Programul calculeaz i vectorii proprii y. Diferena maxim n Ay y este 5.766E-08.
-
a.ch. Noiembrie 2008
10
Programul rulat cu deplasare strategia 3, termin n 9 iteraii, cu diferena maxim 1.228E-12. Pentru acest exemplu, strategia 3 de alegere a deplasrii este cea mai bun.
Exemplu 1 bis: Rafinare Considerm matricea i parametrii din Exemplu-1, i facem rafinarea sistemului propriu cu ns = 4. Se obin rezultatele (valorile proprii sunt listate n format 1pg24.16): Iterations Test Value Tolerance
1 4 4.4408920985E-16 1.776E-15
2 4 2.7194799110E-16 4.441E-16
3 3 0.000000000 8.882E-16
4 4 1.5324295440E-14(stationary)
Eigenvalues
1 14.32950642539381
2 -3.415090280621964
3 4.456959098788065
4 -0.3713752435599113
Diferena maxim n Ay y este acum, 1.110E-15.
Exemplu 2: Matrice nesimetric Considerm matrice:
=
0124051276544321
A
-
a.ch. Noiembrie 2008
11
Programul QR, cu eps = 1E-7, eps_lambda = 1E-8, deplasare strategia 1, produce rezultatele:
Matrix A - Hessenberg form:
1.000000 -5.000000 -1.539516 -1.276672
-6.000000 8.555556 0.1084754 3.698593
0.000000 -5.918166 -2.168879 -1.142756
0.000000 0.000000 -0.1427564 3.613324
Iteration 24:
Reduced Matrix - Block-Triangular Form
4.153893 10.18835 2.745848 2.753036
5.464962 3.096039 3.428772 2.606012
6.2957059E-08 -5.5017566E-09 0.1764519 0.1516238
-2.9451447E-16 2.5737335E-17 1.5892009E-08 3.573617
Eigenvalues
1 -3.855588238007097
2 11.10551971006653
3 0.1764519119325265
4 3.573616616008061
Diferena maxim n Ay y este 4.403E-08. Cu strategia 3, programul termin n 22 iteraii i diferena maxim este 8.755E-08.
Exemplu 2 bis: Rafinare Considerm matricea i parametrii din Exemplu-2, i facem rafinarea sistemului propriu cu ns = 4. Se obin rezultatele:
Iterations Test Value Tolerance
1 4 3.1401849174E-16 4.441E-16
2 4 2.2204460493E-16 1.776E-15
3 5 0.000000000 2.776E-17
-
a.ch. Noiembrie 2008
12
4 5 0.000000000 4.441E-16
Eigenvalues
1 -3.855588220333913
2 11.10551973067810
3 0.1764518729384592
4 3.573616616717359
Diferena maxim n Ay y este 1.776E-15.
2 Problema generalizat
2.1 Problema Problema determinrii valorilor i vectorilor proprii ai unei matrici A, definit de
xAx = , sau de ecuaia
0xIA = )( , (1)
se va numi problema standard.
n Dinamica Structurilor apare urmtoarea problem:
0xMK = )( , (2)
unde: 2 = , iar K i M sunt matrice simetrice i pozitiv definite (K este matricea de rigiditate, M matricea de mas sau de inerie, i este pulsaia proprie).
Problema (2) se zice problema generalizat de valori proprii. Mai general considernd dou matrici nn A i B, unde B este nesingular, problema
generalizat are forma
0xBA = )( (2')
n continuare, vom trata problema (2).
-
a.ch. Noiembrie 2008
13
2.2 Reducerea la problema standard Problema (2) poate fi adus la problema standard (1), prin nmulirea la stnga a ecuaiei (2) cu matricea 1K i notnd MKD 1= , adic,
0xID = )1(
Dezavantajul acestei formulri este c, n general, matricea D nu mai este simetric i astfel, metodele pentru matrici simetrice (de exemplu, Jacobi, iteraii simultane) nu mai pot fi utilizate.
Pe de alt parte:
- Valorile proprii ale unei matrici simetrice sunt bine-condiionate, pe cnd cele ale unei matrici nesimetrice pot fi ru-condiionate.
- La utilizarea metodei QR, numrul de operaii (pe un pas al iteraiei QR) pentru o matrice simetric este mult mai mic dect pentru o matrice general.
innd cont de definirea pozitiv a matricii M, problema generalizat (2) poate fi transformat n problema standard pentru o matrice simetric cum urmeaz.
1. Se face descompunerea Cholesky a lui M:
SSM T= , (3)
unde S este o matrice superior triunghiular. (Descompunerea inferior triunghiular poate fi de asemenea utilizat). Atunci, (2) scris n forma
MxKx = (4)
devine SxSKx T= , i se scrie din nou ca i
)()(1 SxSSxKS T= .
2. Se noteaz
Sxy = (5)
i ecuaia precedent devine
-
a.ch. Noiembrie 2008
14
ySyKS T=1 .
Premultiplicnd cu TTT == SSS )()( 11 (ultima egalitate este o notaie), obinem
IyyKSS = 1T
3. Acum, definim
1= KSSR T (6)
R este o matrice simetric i pozitiv definit v. mai jos. 4. Problema (2) devine problema standard pentru matricea R, i anume,
0yIR = )( (7)
Problemele (7) i (2) au aceleai valori proprii.
5. Dupa rezolvarea lui (7) pentru i y, vectorii proprii ai problemei originale (2) sunt dai (v. (5)) de
ySx 1= (8)
Se verific imediat c matricea R este simetric i pozitiv definit, la fel cum sunt matricile M i K. ntr-adevr, avem:
RKSSSKSKSSR ==== 111)( TTTTTT
Apoi, pentru orice 0x , avem:
0)()()( 111 >=== KuuxSKxSxKSSxRxx TTTTT ,
ntruct xSu 1= este arbitrar i nonzero, la fel ca i x.
Urmeaz c valorile proprii ale lui R sunt reale i pozitive, aa cum trebuie s avem
conform cu 2ii = . Astfel, pentru a rezolva (7), orice metod pentru problema de valori proprii ale unei matrici simetrice i pozitiv definite poate fi aplicat lui R.
Observaie
-
a.ch. Noiembrie 2008
15
Cel mai simplu caz este acela n care matricea M este diagonal (model de mase concentrate), )( imdiag=M . Atunci, avem )( iT mdiag== SS , i
)/1(1 iT mdiag== SS . Astfel, elementele lui R i x sunt date de:
njimm
kr
ji
ijij ,1,, == ,
i
nim
yx
i
ii ,1, == .
Dac M nu este diagonal: notnd 1= SZ , matricea Z se calculeaz uor din
njjj ,1,)()( == eSz prin substituie napoi, ntruct S este superior triunghiular. Apoi, KZZR T= i Zyx = .
Exemplu
Fie matricile:
=
5.200020001
M ;
=
620231
011K
Avem:
=
5.22
1S ;
=
5.2121
11S
Rezult:
=
4.2520525.121
0211R
Vectorii proprii ai lui R sunt calculai analitic, din sistemul (7). Se obine:
-
a.ch. Noiembrie 2008
16
+
=
15.2)1(52
5/2
2 y
Sau, mprind cu prima coordonat
+
=
)15.2(5.2)1(2
1
2 y
Cu acesta, se obine
+
==
15.2)1(2
1
2
1
ySx
Valorile proprii ai lui R se pot calcula rezolvnd ecuaia caracteristic
06.12.69.4)( 23 =++= p . Se obine (n simpl precizie):
=1 3.025604; =2 1.528400; =3 0.3459958.
2.3 Ortogonalitate Se poate arta c vectorii proprii x ai problemei (2), sunt ortogonali relativ la ambele matrici K i M.
1) M-ortogonalitate:
Vectorii y sunt ortogonali, adic,
0)()( =jTi yy , pentru i j .
Substituind Sxy = rezult 0)()( =jTTi SxSx , sau, innd cont de (3),
0)()( =jTi Mxx . (9)
2) K-ortogonalitate:
Vectorii y sunt ortogonali relativ la matricea R, adic
-
a.ch. Noiembrie 2008
17
0)()( =jTi Ryy , pentru i j.
Cum Sxy = , urmeaz c 0)()( =jTTi RSxSx . Din (6) avem KRSS =T , i astfel, obinem
0)()( =jTi Kxx , (10)
care probeaz ortogonalitatea relativ la K. Concluzia (10) se poate obine i din (4), premultipicnd cu Ti)(x i utiliznd (9).
n final, s introducem matricea vectorilor proprii
==
)()2()1(
)(2
)2(2
)1(2
)(1
)2(1
)1(1
)()2()1( ][n
nnn
n
n
n
xxx
xxx
xxx
L
MMMM
L
L
K xxxX,
i s presupunem c X este normalizat, adic 1|||| 2)( =ix . S notm:
iiTi M=)()( Mxx , (11)
iiTi K=)()( Kxx . (12)
Se pot da formule explicite pentru iM i iK . Observai c, n acord cu definirea
pozitiv a lui M i K, avem 0>iM i 0>iK . Mai mult, din (4) obinem ii MK = .
Cu (11) i (12), relaiile de ortogonalitate (9) i (10) pot fi scrise, respectiv, ca
)( iT Mdiag=MXX ,
i
)( iT Kdiag=KXX .
NOTE
-
a.ch. Noiembrie 2008
18
1. Algoritm i Rutine Rutinele pentru calculul pulsaiilor i vectorilor proprii prin metoda de mai sus, se gsesc
n folderele din \ANA\EIGEN\ Biblioteca ANA . Calculul va parcurge urmtorii pai:
1) Calculul matricii R: \Eigen\General_R (sau General_R1) 2) Calculul valorilor i vectorilor proprii ai matricii R:
\Eigen\QR; \Eigen\Jacobi_D; \Eigen\Sim_Iter_D; 3) Regsirea valorilor i vectorilor proprii ai problemei generalizate:
\Eigen\Retrieve_Eigen_from_R Toate rutinele menionate lucreaz n dubl precizie.
2. Note pentru pentru Analiza structurilor - Matricea de rigiditate K va fi furnizat de programul de analiz static cu
ncrcrile date de greutatea maselor. - Matricea de mas M va fi:
Generat direct, pentru mase concentrate; Furnizat de programul de analiz dianmic, pentru mase distribuite. (Acesta va fi rulat pe un interval mic, cu condiii iniiale nenule, i ncrcri nule.)
3. Uniti de msur
Pentru ca pulsaiile s rezulte n 1s trebuie ca elementele matricii M s fie n kg, iar cele
ale matricii K n m
N.
Mai general, n uniti compatibile, astfel ca 2 s rezulte n 2/1 s , iar n s/1 :
Lund fora (F) sau lungimea (L) n alte uniti (dect unitile SI), masa va fi mrime derivat, anume:
2/ tLFM = .
Programul General_R cere, la citirea din fiier a matricilor M i K, introducerea unor
factori care nmulesc aceste matrici. Aceti factori vor asigura transformarea unitilor
folosite la generarea matricilor, n uniti SI (sau, n uniti compatibile).
-
a.ch. Noiembrie 2008
19
4. Pulsaii sau frecvene In particular, General_R permite opiunea de a calcula pulsaii sau frecvene.
Frecvenele sunt legate de pulsaii prin
)2/( pi iif = [rad/s] Perioadele de vibraie sunt:
ii fT
1= [s]