Manipulace s tridiagonálními maticemi
V 1D úloze jsou uzly přirozeně podél zkoumané úsečky – každý má jediného souseda vpravo a jediného
vlevo. Tato vlastnost (neopakující se ve vyšších dimenzích) umožňuje očíslovat uzly (přirozeně) tak, že
matice soustavy bude tridiagonální. Se vzniklou tridiagonální maticí
A =




a1 c1 0 0 . . . 0
b2 a2 c2 0 . . . 0
...
0 . . . 0 0 bN aN




se dá jako s jedním z mála typů laborovat analyticky, například lze nezávisle na její velikosti pomocí
zobecněného kontinuantu
K0 = 1 K1 = a1 Kn = anKn−1 − bn−1cn−1Kn−2
explicitně (rekurzí) vyjádřit její determinant jako det A = KN .
Možnost spočítat v případě potřeby determinant libovolně velké matice rekurzivně (a bez rozvoje
podél řádků/sloupců) představuje obrovské urychlení výpočtu.
Sám determinant tridiagonální matice nestačí pro řešení soustavy MKD, neboť dosazením pravé strany do
jednotlivých sloupců bude tridiagonalita porušena. Přesto existuje způsob, založený na Gaussově eliminaci
(která je obtížnosti O(n3
)), který tridiagonální soustavu umí rozřešit O(n) - TDMA (Tridiagonal
Matrix Algorithm). Tridiagonální soustavu
a1x1 + c1x2 = d1
bixi−1 + aixi + cixi+1 = di i = 2, 3, . . . N − 1
bN xN−1 + aN xN = dN
můžeme převést na tvar
xi + c′
ixi+1 = d′
i i = 1, . . . N − 1
xN = d′
N
kde
c′
1 =
c1
a1
c′
i =
ci
ai − ci−1bi
, i = 2, . . . N − 1
d′
1 =
d1
a1
d′
i =
di − d′
i−1bi
ai − c′
i−1bi
, i = 2, . . . N
Vzniklou soustavu již můžeme zpětným dosazováním snadno rekurzivně vyřešit:
xN = d′
N , i = N − 1, N − 2, . . . 1 : x[i] = d′
i − c′
ixi+1.
(podobně lze efektivně řešit i soustavy, které jsou jen blokově tridiagonální)
Ve speciálním případě digonálně konstantní matice ai = a, bi = b, ci = c lze psát explicitně i vlastní
hodnoty této matice jako
λi = a + 2
√
bc cos
iπ
N + 1
,
explicitně lze v tomto případě nalézt také vlastní vektory.
1
Vlastní hodnoty ve více dimenzích
Při přechodu k vícerozměrným úlohám již není možné uzly číslovat vhodně organizovaným způsobem,
zejména ne v případě MKP. Získané matice tak již nebudou tridiagonální a pro nalezení jejich vlastních
hodnot musíme použít obecnější algoritmy.
V obecném případě pro nalezení vlastních hodnot energie z H − EI použijeme iterativní QR algoritmus:
(QT
Q = I je ortogonální matice, R je horní trojúhelníková matice)
1) nechť H0 ≡ H
2) nalezneme QR dekompozici H0 = Q0R0.
3) určíme Hk+1 = RkQk.
4) loop 2) s Hk+1
Dostatečným počtem kroků postup iteruje k trojúhelníkové matici podobné s výchozí, takže hodnoty z
diagonály jsou přímo vlastní hodnoty řešeného problému.
Vzhledem ke své povaze může být QR rozkald použít také přímo pro solver lineárních soustav.
QR dekompozice s využitím Householderovy transformace
Uvažujme matici H = (hij). QR dekompozice s využitím Householderovy transformace spočívá v následujícím
schematu:
Householderova transformace v roli optimalizace
(získáme horní Hessenbergovu matici, obtížnost (10/3)n3
+ O(n2
))
QR dekompozice
(rozklad na ortogonální a horní trojúhelníkovou matice, obtížnost 6n2
+ O(n)
pro horní Hessenbergerovu matici)
Celkový algoritmus:
1) vezměme první sloupec H, h1j s normou |h1j|
2) vypočteme vektor uj = h1j − |h1j|e1,
kde e1 = (1, 0, 0, . . .)
3) spočteme matici Q[1]ij = Iij − 2
|uj |2 uiuj
Potom H[1] = Q[1]H má v prvním sloupci pod hlavní diagonálou nuly a celý proces můžeme opakovat
iterativně s tím, že se v každém dalším kole zaměříme pouze na zbývající najvětší minor.
Přitom pro násobení Q[k]H[k − 1] musíme Q[k] vždy doplnit shora příslušnou částí jednotkové matice.
Na závěr tvoří Q = Q[1]Q[2] . . . Q[K] a R = QT
H hledanou QR dekompozici H = QR.
V případě Hermiteovské matice jsou vznikající Hessenbergovy matice tridiagonální, a celý algoritmus
Hauseholderovy transformace lze dále optimalizovat. Samotnou otázku tvoří problematika konvergence
(pomalá, jsou-li vlastní hodnoty blízko); řeší se zaváděním posunů spektra.
2
Ukázka kódu pro výpočet QR algoritmu (perl)
for ($cyklus=1;$cyklus<=100;$cyklus++){$kolo=1;
for ($i=$kolo;$i<=$dim;$i++){for ($j=$kolo;$j<=$dim;$j++){$fullQ[$i][$j]=0}};
for ($i=$kolo;$i<=$dim;$i++){$fullQ[$i][$i]+=1};
for ($kolo=1;$kolo<=($dim-1);$kolo++){
$norma=0;
for ($j=$kolo;$j<=$dim;$j++){$u[$j]=$matice[$kolo][$j];$norma+=$u[$j]*$u[$j]};
$u[$kolo]-=sqrt($norma);
$knorma=0;
for ($j=$kolo;$j<=$dim;$j++){$knorma+=$u[$j]*$u[$j]};
for ($i=$kolo;$i<=$dim;$i++){for ($j=$kolo;$j<=$dim;$j++){
$Q[$i][$j]=-2*$u[$i]*$u[$j]/$knorma}};
for ($i=$kolo;$i<=$dim;$i++){$Q[$i][$i]+=1};
for ($i=$kolo;$i<=$dim;$i++){for ($j=1;$j<=$dim;$j++){$nfullQ[$i][$j]=0;
for ($k=$kolo;$k<=$dim;$k++){$nfullQ[$i][$j]+=$fullQ[$k][$j]*$Q[$k][$i]}}};
for ($i=$kolo;$i<=$dim;$i++){for ($j=$kolo;$j<=$dim;$j++){$nmatice[$i][$j]=0;
for ($k=$kolo;$k<=$dim;$k++){$nmatice[$i][$j]+=$Q[$k][$j]*$matice[$i][$k]}}};
for ($j=1;$j<=$dim;$j++){for ($i=1;$i<=$dim;$i++){
$matice[$i][$j]=$nmatice[$i][$j];$fullQ[$i][$j]=$nfullQ[$i][$j]}};
}
print "\n";
for ($i=1;$i<=$dim;$i++){for ($j=1;$j<=$dim;$j++){$puvodni[$i][$j]=0;
for ($k=1;$k<=$dim;$k++){$puvodni[$i][$j]+=$nmatice[$k][$j]*$fullQ[$i][$k]};
$matice[$i][$j]=$puvodni[$i][$j]}};
for ($j=1;$j<=$dim;$j++){for ($i=1;$i<=$dim;$i++){
printf "%1.6f ", $puvodni[$i][$j] };print "\n";}
}
3