2. Matice, soustavy lineárních rovnic Tento učební text byl podpořen z Operačního programu Praha - Adaptabilita
Irena Sýkorová
Některé vlastnosti matic Uvažujme čtvercovou matici A = (aij )n×n . Matice A se nazývá symetrická, jestliže platí A = A′ , tj. aij = aji
pro i, j = 1, . . . , n.
Matice A se nazývá pozitivně definitní, jestliže všechny její hlavní minory jsou kladné, tj.
D1 > 0, . . . , Dn > 0,
kde
¯ ¯ a11 ¯ ¯a Di = ¯ 21 ¯... ¯ ai1
a12 a22 ... ai2
... ... ... ...
¯ a1i ¯ ¯ a2i ¯ ¯ ... ¯ ¯ aii
Matice A se nazývá diagonálně dominantní, jestliže absolutní hodnota prvku na diagonále je větší nebo rovna součtu absolutních hodnot ostatních prvků – buď pro všechny řádky nebo pro všechny sloupce, tj. X |aii | ≥ |aij | pro i = 1, . . . , n. i6=j
Matice je ostře diagonálně dominantní, jsou-li nerovnosti ostré. Příklad 1: Rozhodneme, zda matice A je symetrická, pozitivně definitní, diagonálně dominantní 4 −1 0 A = 2 3 −2 . 0 −2 3 Řešení: Matice A není symetrická, protože a12 6= a21 . K ověření definitnosti vypočítáme příslušné minory ¯ ¯ ¯ ¯ ¯ 4 −1 0 ¯¯ ¯ 4 −1 ¯ ¯ ¯ = 14 (> 0), D1 = 4 (> 0), D2 = ¯¯ D3 = ¯¯ 2 3 −2 ¯¯ = 26 (> 0). 2 3¯ ¯ 0 −2 3¯ Všechny minory jsou kladné, matice provedeme nejprve pro řádky: | | |
A je pozitivně definitní. Podmínku pro diagonální dominanci 4| 3| 3|
> | − 1| < |2| > | − 2|
+ |0| + | − 2| + |0|
ve druhém řádku podmínka neplatí, zkusíme tedy sloupce: |4| > |3| ≥ |3| >
|2| |−1 | | − 2|
+ |0| + | − 2| + |0|
Podmínka je splněná pro všechny sloupce, matice A je diagonálně dominantní, není však ostře diagonálně dominantí (ve druhém sloupci je neostrá nerovnost). MATLAB ověření, zda matice A je symetrická, tedy A = A′ , neboli A − A′ = O (nulová matice): >> A-A’ nebo porovnáním – výsledek: 1 (A je symetrická), 0 (A není symetrická) >> isequal(A,A’) ověření, zda matice A je pozitivně definitní – všechny minory jsou kladné (pro malé matice)
2
>> >> >> >> >>
A(1,1) det(A(1:2,1:2)) det(A(1:3,1:3)) ... det(A)
Vlastní čísla matice Uvažujme čtvercovou matici A = (aij )n×n . Jestliže pro číslo λ (obecně komplexní) a nenulový vektor x ¯ platí A¯ x = λ¯ x, (1) číslo λ se nazývá vlastní (charakteristické) číslo matice A a vektor x ¯ se nazývá vlastní (charakteristický) vektor matice A příslušný vlastnímu číslu λ. Rovnici (1) můžeme upravit, J je jednotková matice řádu n A¯ x = λ¯ x, A¯ x − λ¯ x = o¯,
(A − λJ) x ¯ = o¯.
(2)
Toto je maticový zápis homogenní soustavy n lineárních rovnic o n neznámých s maticí soustavy A − λJ. Hledáme nenulový vektor x ¯, tj. nenulové řešení soustavy, a to existuje, právě když je matice soustavy singulární, tj. det(A − λJ) = 0. (3) Matice A − λJ se nazývá charakteristická matice, polynom p(λ) = det(A − λJ) se nazývá charakteristický polynom (stupně n) a rovnice (3) je charakteristická rovnice. Hledáme kořeny polynomu n-tého stupně a tedy existuje n (reálných nebo komplexních, případně i násobných) vlastních čísel matice A. Spektrální poloměr matice A je největší z absolutních hodnot vlastních čísel ρ(A) = max{|λi |, kde λi je vlastní číslo matice A}.
P o z n á m k a : Známe-li vlastní číslo λ, pak příslušný vlastní vektor je každé nenulové rešení soustavy (2). Příklad 2: Vypočítáme vlastní čísla a spektrální poloměr matice ¶ µ 1 2 . A= 3 2 Řešení: Nejprve určíme chrakteristickou matici µ ¶ µ 1 2 λ A − λJ = − 3 2 0 pak charakteristický polynom, neboli ¯ ¯1 − λ det(A − λJ) = ¯¯ 3
0 λ
¶
=
µ
1−λ 2 3 2−λ
¶
,
¯ 2 ¯¯ = (1 − λ)(2 − λ) − 6 = λ2 − 3λ − 4. 2 − λ¯
Vlastní čísla matice jsou kořeny charakteristického polynomu, tj. hledáme řešení charakteristické rovnice (kvadratické) λ2 − 3λ − 4 = 0, tedy charakteristická čísla matice A jsou λ1 = 4, λ2 = −1.
3
Spekrální poloměr je největší z absolutních hodnot vlastních čísel, tj. max{|λ1 |, |λ2 |} = max{| 4 |, | − 1|} = max{4, 1} = 4. Spekrální poloměr je ρ(A) = 4. Platí: a) Je-li matice symetrická, všechna její vlastní čísla jsou reálná. b) Matice je symetrická a pozitivně definitní, právě když všechna její vlastní čísla jsou kladná. MATLAB vlastní čísla matice A: >> eig(A) spektrální poloměr matice A (největší z absolutních hodnot vlastních čísel): >> max(abs(eig(A)))
Norma matice Podobným způsobem jako u vektorů můžeme definovat i normu matice. Norma matice je takové zobrazení Rn×n do R, které splňuje tyto podmínky: (i) (ii) (iii) (iv)
kAk > 0 pro každou nenulovou čtvercovou matici a kAk = 0 pouze pro A = O, kc · Ak = |c| · kAk pro každé reálné číslo c, kA + Bk ≤ kAk + kBk pro každé dvě čtvercové matice stejného řádu. kA · Bk ≤ kAk · kBk pro každé dvě čtvercové matice stejného řádu.
Je důležitá souvislost mezi normou vektoru a normou matice. Řekneme, že norma vektoru a norma matice je konzistentní s normou vektoru, jestliže pro každou matici A a vektor x ¯ platí kA¯ xk ≤ kAk · k¯ xk. Normu matice můžeme také definovat pomocí normy vektoru. Je-li dána norma vektoru k · k, pak platí kA¯ xk = sup kA¯ xk (4) kAk = sup k¯ x k x ¯6=o ¯ k¯ xk=1 a říkáme, že tato norma matice je generovaná normou vektoru k · k. Norma matice generovaná normou vektoru je s touto normou konzistentní. Nejčastěji užívané normy matice jsou: Pn a) kAk∞ = max j=1 |aij |, tzv. řádková norma, 1≤i≤n Pn b) kAk1 = max i=1 |aij |, tzv. sloupcová norma, 1≤j≤n ³P ´1 n Pn 2 2 c) kAkF = , tzv. Frobeniova norma, někdy též Schurova, j=1 aij i=1 p ′ d) kAk2 = ρ(A A), tzv. spektrální norma, někdy též euklidovská.
4
Příklad 3: Vypočítáme všechny normy matice A =
µ
¶ 1 −5 . 2 6
Řešení: a) Při počítání řádkové normy sečteme absolutní hodnoty prvků v každém řádku a z těchto součtů vybereme maximální kAk∞ = max{|1| + | − 5|, |2| + |6|} = max{6, 8} = 8. b) Pro určení sloupcové normy sečteme absolutní hodnoty prvků v každém sloupci a z těchto součtů vybereme maximální kAk1 = max{|1| + |2|, | − 5| + |6|} = max{3, 11} = 11. c) Frobeniova norma je druhá odmocnina součtu druhých mocnin všech prvků matice p √ √ kAkF = (1)2 + (2)2 + (−5)2 + (6)2 = 1 + 4 + 25 + 36 = 66 = 8, 124.
d) K určení spektrální normy potřebujeme nalézt největší z absolutních hodnot vlastních čísel matice A′ A ¶ ¶ µ ¶ µ µ 5 7 1 −5 1 2 ′ , = · AA= 7 61 2 6 −5 6
pak kořeny chakteristického polynomu jsou (5 − λ)(61 − λ) − 49 = 0
⇔
λ2 − 66λ + 256 = 0
⇔
λ1 = 61, 8617, λ2 = 4, 1383.
Spektrální norma je odmocnina ze spektrálního poloměru matice A′ A p ρ(A′ A) = 61, 8617 ⇒ kAk2 = 61, 8617 = 7, 8652.
Výpočet spektrální normy pro větší matice je poměrně pracný. MATLAB – normy matice A: řádková: >> norm(A,inf) >> max(sum(abs(A’))) sloupcová: >> norm(A,1) >> max(sum(abs(A))) Frobeniova: >> norm(A,’fro’) >> sqrt(sum(diag(A’* A))) spektrální: >> norm(A,2) >> norm(A)
P o z n á m k a : Frobeniova norma matice není√generovaná euklidovskou normou vektoru, protože pro jednotkovou matici J řádu n platí kJkF = n, ale podle (4) je kJkF = 1. Euklidovskou normou vektoru je generovaná spektrální norma matice. Pro každou normu matice, která je konzistentní s nějakou normou vektoru, platí ρ(A) ≤ kAk.
5
Podmíněnost matic Dříve než se budeme věnovat studiu metod řešení soustav lineárních rovnic, je potřeba zmínit tzv. podmíněnost matice. V podstatě jde o to, jak „citliváÿ je matice soustavy vzhledem k chybám ve vstupních datech i k zaokrouhlovacím chybám v průběhu výpočtu. Příklad 4: Uvažujme soustavu lineárních rovnic 2x + 6y = 8 2x + 6, 0001y = 8, 0001, jejímž řešením je x = 1, y = 1. Pokud provedeme v koeficientech malou změnu (řádově 104 ), dostaneme soustavu 2x + 6y = 8 2x + 5, 9999y = 8, 0002 s řešením je x = 10, y = −2. Tedy malé změny koeficientů matice a složek vektoru pravých stran způsobí velké změny (řádově jednotky) v řešení. Matice se nazývá dobře podmíněná, jestliže relativně malé změny v koeficientech způsobí relativně malé změny v řešení. Matice se nazývá špatně podmíněná, jestliže relativně malé změny v koeficientech způsobí relativně velké změny v řešení. Uvažujeme nyní soustavu lineárních rovnic A¯ x = ¯b
(5)
s regulární maticí soustavy A. Označíme-li x ¯∗ přesné (teoretické) řešení soustavy A¯ x = ¯b a x ¯c je x∗ k c −¯ přesné řešení porušené soustavy (A + δA)¯ x = ¯b + δ¯b, pak odhad relativní chyby řešení k¯xk¯ závisí x∗ k přímo na součinu K(A) = kAk · kA−1 k, kde norma matice je generovaná normou vektoru. Číslo K(A) se nazývá číslo podmíněnosti matice A. Čím je větší, tím je matice hůř podmíněná a tím je větší odhad relativní chyby řešení. µ ¶ 2 6 P o z n á m k a : Matice A = z předchozího příkladu má číslo podmíněnosti velké 2 6, 0001 K(A) = 4 · 105 . MATLAB číslo podmíněnosti matice A: >> cond(A)
Přímé metody řešení soustav lineárních rovnic Nyní se budeme věnovat numerickému řešení soustavy lineárních rovnic A¯ x = ¯b, kde matice soustavy A = (aij ) je reálná regulární matice řádu n, ¯b = (b1 , b2 , . . . , bn )′ je sloupcový vektor pravých stran.
6
Metody řešení soustav lineárních rovnic můžeme rozdělit do dvou skupin: metody přímé a metody iterační. Pomocí přímých metod dostaneme po konečném počtu kroků přesné řešení soustavy, metodami iteračními získáme posloupnost vektorů, která konverguje k přesnému řešení. Ve skutečnosti ale vždy dostaneme pouze určitou aproximaci řešení; u přímých metod je to způsobeno zaokrouhlovacími chybami, u metod iteračních tím, že můžeme provést vždy jen konečný počet kroků. Volba metody závisí na konkrétní soustavě, musíme vzít v úvahu, zda je matice soustavy malá nebo velká, zda obsahuje hodně nulových prvků, tzn. je řídká nebo zda má nějaké speciální vlastnosti – např. je třídiagonální. Mezi přímé metody patří Gaussova eliminační metoda, Jordanova eliminační metoda, metoda využívající inverzní matici. Známe-li inverzní matici A−1 , můžeme počítat řešení ze vztahu x ¯ = A−1¯b, ale vlastní výpočet matice A−1 není příliš výhodný. Gaussova eliminační metoda spočívá v tom, že nejprve rozšířenou matici soustavy (A|b) převedeme pomocí ekvivalentních úprav do tvaru (U |y), kde U je horní trojúhelníková matice – tzv. přímý chod. Místo původní soustavy pak řešíme soustavu Ux ¯ = y¯, ze které se zdola snadno dopočítají složky neznámého vektoru x ¯ – tzv. zpětný chod. V přímém chodu provádíme v (n−1) krocích nulování prvků pod diagonálou v 1. až (n−1). sloupci. Používáme pouze takovou úpravu, že v k-tém kroku násobíme k-tý řádek vhodnou konstantou a tento násobek přičteme k ostatním řádkům. Označíme-li A(0) = A, b(0) = b, počítáním podle vzorce (pro k=1, . . . , n-1, i=k+1, . . . , n, j=k, . . . , n) (k) aij
=
(k−1) aij
−
(k−1) (k−1) akj aik (k−1) akk
(k−1)
s podmínkou akk 6= 0. Je-li příslušný diagonální prvek nulový, musíme tento řádek vyměnit s některým ze zbývejících, tak aby na diagonále bylo nenulové číslo. Případně, kvůli omezení šíření zaokrouhlovacích chyb, je možné provést pivotaci – výběr prvku v absolutní hodnotě největšího v daném sloupci, který výměnou řádků přemístíme na diagonální pozici. Ve zpětném chodu se neznámé počítají odzadu podle vzorců yn xn = unn n X 1 uij xj pro i = n − 1, . . . , 1. yi − xi = uii j=i+1 Příklad 5: Gausovou eliminační metodou vyřešíme soustavu x1 + 2x2 + 2x3 =
2
2x1 + 2x2 + 3x3 = 5 −x1 − x2 = −1.
(6)
Řešení: Koeficienty a pravé strany zapíšeme do rozšířené matice soustavy, kterou převedeme do horního trojúhelníkového tvaru. Provádíme takové úpravy, že v k-tém kroku násobíme vždy k-tý řádek vhodnou konstantou a tento násobek přičteme k ostatním řádkům tak, abychom vynulovali prvky pod diagonálou (konstanty jsou zapsané vpravo vedle příslušného řádku). 1 2 2 1 2 2 1 2 2 2 2 2 (−2) (1) 2 2 3 ∼ 0 −2 −1 5 1 ( 12 ) ∼ 0 −2 −1 1 3 3 −1 −1 0 −1 0 1 2 1 0 0 2 2
7
Odtud snadno dopočítáme x3 = 1, x2 = −1, x1 = 2.
Podobná je Jordanova eliminační metoda, kde se matice soustavy A převádí ma matici jednotkovou, řešení pak je ve sloupci pravých stran. Z Gaussovy eliminace můžeme odvodit i tzv. LU rozklad matice. Dříve než popíšeme princip LU rozkladu, ukážeme, že některé elementární úpravy matice lze provádět násobením dané matice zleva speciálními regulárními maticemi. P o z n á m k a : Jestliže v dané matici M chceme m násobek i-tého řádku přičíst k j-tému řádku, můžeme tuto úpravu vyjádřit jako součin matic; matici M nádobíme zleva maticí V , která vznikne z jednotkové matice přidáním konstanty m na pozici vji . Pokud je j > i, je matice V dolní trojúhelníková. Příklad 6: Matici M upravíme tak, že čtyřnásobek prvního řádku přičteme ke třetímu (i = 1, j = 3, m = 4). −1 2 0 1 0 0 2 2 V = 0 1 0 M = 3 4 −5 2 4 0 1 −1 2 0 −1 2 0 1 0 0 2 2 = 3 2 2 V · M = 0 1 0 · 3 0 3 2 4 −5 2 4 0 1 P o z n á m k a : Jestliže v dané matici M chceme vyměnit i-tý a j-tý řádek, můžeme tuto úpravu vyjádřit jako součin matic; matici M nádobíme zleva maticí P , která vznikne z jednotkové matice záměnou i-tého a j-tého řádku. Matice P se nazývá permutační matice. Příklad 6: V matici M vyměníme první a druhý řádek (i = 1, j = 2). −1 2 0 0 1 0 M = 3 2 2 P = 1 0 0 4 −5 2 0 0 1 3 2 2 −1 2 0 0 1 0 2 0 2 2 = −1 P · M = 1 0 0 · 3 4 −5 2 4 −5 2 0 0 1 Budeme-li používat LU rozklad matice, znamená to, že matici soustavy A chceme vyjádřit jako součin A = LU, kde L je dolní trojúhelníková matice s jedničkami na diagonále a U je horní trojúhelníková matice. Soustavu (5) pak řešíme jako dvě soustavy s trojúhelníkovými maticemi L¯ y = ¯b Ux ¯ = y¯. Trojúhelníkové matice L a U můžeme získat postupnými úpravami matice A, což můžeme symbolicky vyjádřit jako (A|J) ∼ . . . ∼ (U |L′ ),
kde L′ je dolní trojúhelníková a platí L′ A = U . Při tomto způsobu výpočtu však smíme použít jen takovou elementární úpravu, kdy přičítáme násobek určitého řádku k řádku, který je pod ním. Ze vztahu L′ A = U snadno odvodíme A = (L′ )−1 U = LU. Provedeným úpravám odpovídá postupné násobení matice A maticemi Vk L′ A = Vk . . . V1 A = U
⇒
A = (L′ )−1 U = (Vk . . . V1 )−1 U = (V1−1 . . . Vk−1 ) U = LU.
Každá matice L má na diagonále jedničky a pod diagonálou příslušné koeficienty z úprav při Gaussově eliminaci, ale s opačnými znaménky.
8
Příklad 7: Vyřešíme soustavu (6) pomocí LU rozkladu. Řešení: Matici soustavy (6) rozložíme
1 0 0 1 1 0 · 0 A = LU = 2 −1 − 21 1 0
2 2 1 2 −2 −1 = 2 2 3 0 −1 −1 2
2 3. 0
Pak hledáme řešení soustavy L¯ y = ¯b 1 0 0 2 y1 2 1 0 y2 = 5 , −1 − 12 1 −1 y3
kterou snadno vyřešíme: z první rovnice je y1 = 2, dosazením do druhé dostaneme y2 = 1, ze třetí y3 = 23 , tedy 2 y¯ = 1 . 3 2
Nyní zbývá vyřešit soustavu U x ¯ = y¯, což je zpětný chod Gaussovy eliminace.
1 2 2 2 x1 0 −2 −1 x2 = 1 3 3 0 0 x3 2 2
⇒
2 x1 x2 = −1 1 x3
I když matice A je regulární, může se stát, že se během úprav objeví nulový diagonální prvek. Pak je potřeba vyměnit řádky, což v tomto případě není povolená úprava. Je tedy třeba prohodit řádky původní matice A tak, aby se nová matice P A dala se rozložit na součin LU , tj. P A = LU. Řešení původní soustavy pak probíhá takto: L¯ y = P ¯b Ux ¯ = y¯. Výhoda LU rozkladu se projeví, pokud řešíme více soustav se stejnou maticí. Pak se nejpracnější část výpočtu, tj. výpočet matic L a U , provádí pouze jednou. MATLAB – přímé metody řešení soustav lineárních rovnic (A – matice soustavy, b – sloupcový vektor pravých stran): příkazem Matlabu: >> x = A \ b pomocí inverzní matice: >> A1 = inv(A) >> x = A1 * b pomocí LU rozkladu (řešíme Ly = P b, U x = y, kde P je permutační matice): >> [L, U, P] = lu(A) >> x = U \ ( L \ (P*b))
9
Dodatky Je-li matice A symetrická a pozitivně definitní, pak existuje horní trojúhelníková metice U tak, že platí A = U ′ U. (7) Tomuto vyjádření se říká Choleského rozklad. Vztah (7) bývá někdy zapsán také ve tvaru A = LL′ , kde L je dolní trojúhelníková matice. MATLAB Choleského rozklad matice A (funkce vrací horní trojúhelníkovou matici): >> chol(A) řešení soustav lineárních rovnic : >> x = U \ ( U’ \ b) V některých případech potřebujeme vyřešit soustavu s maticí ve speciálním tvaru, maticí třídiagonální, tj. maticí ve tvaru a11 a21 0 . . . 0 0
0
a12 a22 a32 .. .
0 a32 a33 .. .
0 0 a34 .. .
0 0 0 .. .
... ... ... .. .
... ... ... .. .
0 0 0 .. .
... ... ...
... ... ...
0 0 0
an−2,n−3 0 0
an−2,n−2 an−1,n−2 0
an−2,n−1 an−1,n−1 an,n−1
0 an−1,n an,n
.
Taková matice je řídká, má málo nenulových prvků a hodně nulových (pro velká n). Při numerickém řešení se používá speciální postup, který umožňuje uložit do paměti pouze nenulové prvky.
10
