Řešení soustav lineárních rovnic

Řešení soustav lineárních rovnic

Obsah

Řešení soustav lineárních rovnic přímé metody Gaussova eliminace LU rozklad iterační metody prostá Jacobiova Gaussova-Seidelova gradientní metody největšího spádu sdružených gradientů

2

Řešení soustav lineárních rovnic • Úloha. Dána soustava rovnic ve tvaru Ax = b. Určete neznámý vektor x. • Řešitelnost. Je li matice A regulární (det A různý od nuly), existuje jediné řešení. (podrobněji viz Skripta Matematika 1, Frobeniova věta)

3

Přímé a iterační metody Přímé metody

• Za konečný počet kroků určíme přesné řešení (neuvažujeme li zaokrouhlovací chyby výpočtu) • Příklady: – Gaussova eliminační metoda – LU rozklad matice soustavy

Iterační metody • Volíme počáteční aproximaci - vektor x(0) , a postupně vypočítáváme vektory x(1), x(2), … , které se od přesného řešení x* liší tím méně, čím více iteračních kroků provedeme • Příklady: • Maticové iterační metody – Prostá, Jacobiova, GaussovaSeidelova

• Gradientní metody – největšího spádu – sdružených gradientů 4

Řešení soustav lineárních rovnic přímé metody Gaussova eliminační metoda. Ekvivalentními úpravami (násobení řádu nenulovým číslem a přičítání k řádku násobku jiného řádku) převedeme soustavu rovnic na ekvivalentní soustavu s trojúhelníkovou maticí soustavy (přímý chod). Z poslední rovnice (ve které je pouze jedna neznámá) vypočítáme poslední neznámou, kterou dosadíme do zbývajících rovnic, tak pokračujeme k první rovnici (zpětný chod). 5

Gaussova eliminační metoda - přímý chod • •

Označíme A(0) (rozšířenou) matici soustavy - A (nebo A|b) Krok 1: vytvoříme matici A(1) – – –

•

Krok 2: vytvoříme matici A(2) – – –

•

•

Opíšeme první, druhý, …, i-tý řádek matice A(i-1), Ke (i+1) ,…, poslednímu řádku matice A(i-1) přičteme takové násobky i-tého řádku matice A(i-1), abychom v i-tém sloupci získali nulové prvky tj. k (i+1)-mu řádku přičítáme (-a i+1,i)/a ii násobky i.řádku, …, k poslednímu (n-tému) řádku přičítáme (-ani)/a ii násobky i.řádku.

Krok poslední (n-1): vytvoříme matici A(n-1) – –

•

Opíšeme první a druhý řádek matice A(1), Ke třetímu ,…, poslednímu řádku přičteme takové násobky druhého řádku matice A(1), abychom ve druhém sloupci získali nulové prvky tj. ke třetímu řádku přičítáme (-a32)/a22 násobky 2.řádku, ke čtvrtému řádku přičítáme (-a42)/a22 násobky 2.řádku,…, k poslednímu (n-tému) řádku přičítáme (-an2)/a22 násobky 2.řádku. (V LU rozkladu(viz dále) budou právě tyto násobky s opačným znaménkem tvořit 2.sloupec matice L)

Krok i: vytvoříme matici A(i) – –

•

opíšeme první řádek matice A(0), ke druhému, …, poslednímu řádku matice A(0) přičteme takové násobky prvního řádku, abychom v prvním sloupci získali nulové prvky, tj. ke druhému řádku přičítáme (-a21)/a11 násobky 1.řádku, ke třetímu řádku přičítáme (-a31)/a11 násobky 1.řádku,…, k poslednímu (n-tému) řádku přičítáme (-an1)/a11 násobky 1.řádku. (V LU rozkladu(viz dále) budou právě tyto násobky s opačným znaménkem tvořit 1.sloupec matice L)

Opíšeme první, druhý, …, (n-1) řádek matice A(n-2), K poslednímu (n-tému) řádku přičteme takové násobky (n-1) řádku, abychom v (n-1) sloupci získaly nulové prvky.

Poznámka 1: stane li se, že v některém kroku (k) bude prvek a kk nulový, vyměníme řádek k s řádkem m, m>k. Podmínka , že matice A je regulární zaručuje, že je to možné. Poznámka 2: aby vliv zaokrouhlovacích chyb byl menší, je možné v každém kroku i nejprve vyměnit řádek i s takovým řádkem m (m > k), ve kterém je prvek a mi největší v absolutní hodnotě. 6

Gaussova eliminační metoda - zpětný chod •

Matice A(n-1), vytvořená přímým chodem vypadá: a a (0) a (0)  a (0) b ( 0 )  (horní index v závorce je pro připomenutí,  (1)  na kterém kroku byl tento prvek vypočten) (1) (1) (1) a a a b  0 11

 0(1)   0(1)  (1) 0

12

13

1n

22

23

2n

( 2)

0 0( 2)

 

0( 2)



 a n1 ,n1 0( n 1)

a (nn1,n2 ) ann

1

  bn( n1 2 )   bn( n 1)  2

Poslednímu řádku odpovídá tvar poslední rovnice : ann xn = bn • Z poslední rovnice vypočteme poslední neznámou xn = bn / ann, kterou dosadíme do (všech) rovnic n-1…1 • Z předposlední rovnice, která po dosazení vypadá: an-1,n-1 xn-1+ an-1,n xn = bn-1 vypočítáme předposlední neznámou xn-1=(bn-1- an-1,nxn)/an-1,n-1 …

• • •

Ze třetí rovnice vypočteme x3 = (b3- a3,nxn- a3,n-1xn-1 - … -a3,4x4)/a3,3 Z druhé rovnice vypočteme x2 = (b2- a2,nxn- a2,n-1xn-1 - … -a2,3x3)/a2,2 Z první rovnice vypočteme x1 = (b1- a1,nxn- a1,n-1xn-1 - … -a1,2x2)/a1,1

7

Gaussova eliminační metoda: soustava s třídiagonální maticí • Uvažujme soustavu rovnic:

0 0 0   x1   b1   a11 a12 0       a a a 0 0 0  21 22   x2   b2  23  0 a a33 a34 0 0   x3   b3  31       0 0 a a a 0    x4   b4  41 44 45  0   x  b  0 0 a a a 51 55 56    5  5  0 0 0 0 a61 a66   x6   b6  

Potom v každém kroku (i) přímého chodu přičítáme násobky řádku i pouze k jednomu řádku i+1. Matice A(0)… A(n-1) můžeme zapisovat jako trojice vektorů u=(u1, u2,…,un-1) = (a12,a23,…,an-1,n), d=(d1,d2,…,dn) = (a11,a22,…,an,n), l=(l 1, l 2,…, l n-1) = (a21,a32,…,an,n-1) a vektoru pravé strany b=(b1,b2,…,bn) Úpravy přímého chodu (které mění pouze vektory d a b) lze zapsat: di=di+ui-1(- l i-1/di-1), bi=bi+bi-1(- l i-1/di-1),i=2,…,n Zpětný chod je vyjádřen: xn=bn/dn, xi= (bi-ui xi+1) /di,i=n-1,…,1

8

• Příklad 1.

Gaussova eliminační metoda příklady Příklad 2.  4 1 0 2 0  0 0  0 

 6 3 1 x   1        2 2 1 y    3   3 2 1 z   2       Přímý chod : (0)

ozn.

A

Krok 1 :

6 3 1  :2 2 1 3 2 1 

2.ř  1.ř * ( 62 ) 3.ř  1.ř * ( 63 )

1  3 2  (1)

A

Krok 2 : 3.ř  2.ř * ( 21 ) A

6 3 1  :  0 1 23 0 0 1 6 

1  8 3 3 2 1  8 3 1 6

Zp.chod : posl. rovnice : 16 z  16  z  1, druhá rovnice : první rovnice :

1 4 2 0 0

0 0 0 1 4 2

0  x1   6      0  x2   13  0  x3   20       0  x4   27    1  x5   34  4  x6   34 

Označíme vektory

6 3 1  :  0 1 23 0 1 1 2 2 

( 2)

4 2 0 0 0

0 0 1 4 2 0

y  23 z  83  y  2, 6 x  3 y  z  1  x  1

 4 6      1  2      4  13  1 2          4       20  u  1 ,d  ,l  2 ,b       4  27  1    2  4  34   1          2  4  34      Přímý chod : Zpětnýchod :  4   6  1       7 10  2     2  3   24 7    1007  , b    x  d  56 41  12   3   4  140   946  5 41 41        239   1434  6  70   70   

9

Gaussova eliminační metoda a LU rozklad matice •

•

V situaci, kdy je potřeba řešit několik (mnoho) soustav rovnic se stejnou maticí soustavy, ale různými pravými stranami, je vhodné jednou rozložit matici soustavy A na součin trojúhelníkových matic L a U, a potom pro každou pravou stranu řešit 2 soustavy rovnic s trojúhelníkovými maticemi (což je méně pracné): Uy=b (tím aplikujeme úpravy přímého chodu na vektor b) a Lx=y (odpovídá zpětnému chodu) Matice L slouží k „zapamatování“ úprav přímého chodu, U odpovídá matici A(n-1), která vzniká přímým chodem Gaussovy eliminační metody (za předpokladu, že jsme a) v kroku i přičítali násobky i-tého řádku (a ostatní řádky ničím nenásobili) a b) neměnili pořadí řádků matice). Horní index (i) označuje prvek matice, získaný v kroku i přímého chodu Gaussovy eliminace.  1  a (0)  21( 0 )  a11  a (0) L   31( 0 )  a11    a (0)  n1( 0 )  a11

0

0



1

0



1







(1)

a32 (1) a22  (1) an1 (1) a22

( n2)

an1  ( n2) ann

0   a11( 0 ) 0    0  0 , U   0        0 1  

(0)

a12 (1) a22 0  0

( 0) (0) a13  a1n   (1) (1) a23  a2 n  ( 2)  ( 2) a33  a3n     ( n 1)   0  ann 

10

LU rozklad matice - příklad • Příklad 3. Řešte soustavy rovnic Ax=b, Az=c užitím LU rozkladu  6 3 1 1  45  matice A.       A   2 2 1, b   3 , c   24   3 2 1  2  28       

•

a)LU rozklad: v prvním kroku násobíme přičítáme ke druhému řádku (-2/6)násobky prvního řádku a ke třetímu (-3/6) násobky prvního řádku(viz příklad 1), proto v prvním sloupci hledané matice L budou tyto násobky s opačným znaménkem; ve druhém kroku přičítáme ke třetímu řádku (-1/2) násobky druhého řádku, proto ve druhém sloupci matice L bude tento násobek s opačným znaménkem. Hledaná matice U odpovídá upravené matici soustavy. -3/6-

-2/6

6 3 1 6 3 1  6 3 1       2 2  1 3   0 1 3  -1/2  0  2 2 1  0 1 1  0 0 1   3 2 1   2 2 6  

1 0 0 6 3 1 2    2 L   6 1 0 ,U   0 1 3   3 1 1 0 0 1  6 6 2   11

Řešení soustavy užitím LU rozkladu • b) pomocí matice L upravíme vektor pravé strany b (řešíme Ly=b) a určíme neznámý vektor x (řešíme Ux=y )==>x3 =1,x2=2,x1=-1   L y b      1 0 0   y1   1   y1  1  2       6 1 0    y2    3   y2  8 / 3   3 1 1   y   2  y  1/ 6 3 6 2   3  

U y     x1  1  6 3 1   x1   1        2 8 0 1   x2  2   3   x2    3  0 0 1   x   1  x 1  3  6 6  3 

• c)stejným způsobem upravíme vektor c: řešíme Lw=c a určíme neznámý vektor z: řešíme Uz=w ==>z3=6,z2=5,z1=4 

 L c w      1 0 0   y1   45   y1  45  2       6 1 0    y2    24   y2  9   3 1 1   y   28   y  1 3 6 2   3  



U w     z1  4  6 3 1   z1   45         0 1 23    z2    9   z2  5  0 0 1   z   1  z  6  3   6  3 

12

Některé operace s maticemi v MATLABu • • • • • • • • • • • • • •

Malá a velká písmena se rozlišují Proměnné není potřeba deklarovat Operace s maticemi (+,-, *,/,…) jsou DEFINOVÁNY Zadání matice: po řádcích, prvky v řádku odděleny mezerou, řádky odděleny středníkem: A=*a11 a12 a13; … a31 a32 a33+ Vektor řádek (tj. matice 1xn): b=(b1 b2 …bn) Vektor sloupec(tj. matice nx1):c=(c1;c2; …cn) Získat (změnit) prvek matice: z=A(1,2) … A(1,2)=4 Násobení matic(rozměry musí odpovídat!): d=A*c výsledný vektor sloupec, f=A*b nelze, q=b*A výsledný vektor řádek Transponování: b’ tedy f=A*b’ lze, výsledný vektor sloupec Řešení soustavy Ax=b: x=A \ b LU rozklad: [L U]=lu(A) (pozor, algoritmus je jiný---) Určení determinantu matice: d=det(A) Určení vlastních čísel matice: vl_cisla=eig(A) Inverzní matice: Ainv=inv(A) 13

Kontrola výsledků v MATLABu (Gaussova eliminace, LU rozklad) • • • • • • • • •

A=[6 3 1; 2 2 1; 3 2 1] b=[1; 3; 2] x_gauss=A\b [L,U]=lu(A) b_upr=L\b x_LU=U\b_upr b=[45; 24; 28] b_upr=L\b x_LU=U\b_upr

//zadat matici A //zadat vektor b //vypocet x Gauss.el. //vytvoreni LU rozkladu //reseni Lb_upr=b //reseni LU rozkladem //zadani jine prave strany //reseni Lb_upr=b //reseni LU rozkladem 14

Řešení soustav lineárních rovnic maticové iterační metody • Soustavu rovnic Ax=b převedeme na soustavu ve tvaru x=Ux+v • Volíme x(0), počítáme x(1) = Ux(0)+v, x(2) = Ux(1)+v,… x(k+1) = Ux(k)+v

• Způsob, jak vytvoříme matici U a vektor v je určen iterační metodou • Aby posloupnost vektorů x(1),…,x(k) konvergovala k přesnému řešení, musí být splněny podmínky konvergence iterační metody 15

• Pojmy (A : čtvercová matice , prvky – reálná čísla) – ostře diagonálně dominantní matice (ODD) po řádcích : když i = 1,…,n platí |aii|>| aij |, j =1,…,n, j≠i po sloupcích: když i = 1,…,n platí |aii|>| aji |, j =1,…,n, j≠i – A je ODD, když je ODD po řádcích nebo po sloupcích – symetrická matice A=AT tj. aij=aji, i,j = 1,…,n – pozitivně definitní matice jsou li všechny hlavní minory kladné – norma matice • řádková ||A|| = max i  | aij |, j =1,…,n • sloupcová ||A||1 = max j  | aij |, i =1,…,n • Euklidovská ||A||1 =( | aij |2 ) 1/2 , i =1,…,n, j =1,…,n

– vlastní čísla matice – řešení charakteristické rovnice det(A-E)=0

– spektrální poloměr matice je největší z absolutních hodnot vlastních čísel tj. (A) = max i | i | – absolutní hodnota komplexního čísla: z=a+ib, |z| = (a2+b2) 1/2 16

Příklady (vlastnosti matic) •

Dány matice, vektor. Určete normy matic a vektoru, spektrální poloměr matic. Jsou matice symetrické, pozitivně definitní, ostře diagonálně dominantní?

 4 1 0 0.5     0.4  0.1 0   1       0 5 0 0  A B   0.1 0.2 0.4  c   0.5  0 2 2 0   0  0.4 0.4   0.7        0 0 0 1    Normy

4

4

4

4

|a1 j | |a2 j | |a3 j | |a4 j | j 1  j1 j1 j1  || A ||   max{ 5.5 , 5, 4, 1 }  5.5 4

4

4

1.r . 2. r .  3.r .   || B ||   max{0.5, 0.7, 0.8 }  0.8

4

    | ai 4 | i 1 1  i1 i1 i || A ||1  max{ 4, 8, 2, 1.5 }  8

1. sl .  2. sl .  3. sl .  || B ||1  max{0.5, 0.7, 0.8 }  0.8

|| A || 2  4 2  12  ...12  51.25  7.2

|| B || 2  0.4 2  0.12  ...0.4 2  0.7

| ai 1 |

| ai 2 |

| ai 3 |

|| c ||   max{1,0.5,0.7}  1, || c ||1  1  0.5  0.7  2.2, || c || 2  1  0.25  0.49  1.72 17

Příklady (vlastnosti matic) spektrální poloměr det( A  E )  ( 4   )(5   )(2   )(1   )  1  4, 2  5, 3  2, 4  1 spektrální poloměr :  (A) | 2 | 5 det(B  E )  (0.4   )(0.2   )(0.4   )  0.16 (0.4   )  0.01(0.4   ) det(B  E )  (0.4   )(2  0.6  0.25) 0.6  i 0.64  0.3  i  0.4 | 2 || 3 | 0.33  0.4 2  0.5 2 spektrální poloměr :  (B)  0.5

1  0.4 2,3 

symetrie A není symetrická (např.) : a12  a21, B není symetrická : 0.1  0.1,0.4  0.4

pozitivní definitnost A : minory : a11  0,

4 0

1  20  dál nemusíme počítat, není pozitivně definitní 5

B : minory : b11  0,

0.4  0.1  0.09  0, det(B )  0.1  0  je pozitivně definitní 0.1 0.2

ostrá diagonální dominantnost A : po řádcích není (nesplněno ve 3.ř.), po sloupcích je  A je ODD B : po řádcích není (nesplněno ve 2.,3.ř.), po sloupcích není (2.,3.sl.)  B není ODD 18

Kontrola výsledků v MATLABu - vlastnosti matice • • • • •

• • • • • • • • •

•

A=[0.4 -0.1 0; 0.1 0.2 0.4; 0 -0.4 0.4] //zadani matice A (3x3) po radcich vlastni_cisla=eig(A) //vektor vlastnich cisel rho=max(abs(vlastni_cisla)) //spektralni polomer isequal(A,A’) //vysledek: 1 - kdyz A je symetricka, 0 - neni A(1,1)>0 & det([A(1,1) A(1,2); A(2,1) A(2,2)])>0 & det(A)>0 //vysledek: 1 - kdyz A je pozitivne definitni, 0 – neni //ODD D=abs(diag(A)) // D: vektor absolutnich hodnot A(i,i) [M,N]=size(A) //M:pocet radku, N: pocet sloupcu for(i=1:1:M) plati=D(i)>sum(abs(A(i,:)))-D(i); end //plati je vektor z 0 a 1, 1- kdyz podminka splnena ODD=isequal(plati, ones) //vysledek: 1 – ODD po radcich, 0 neni for(i=1:1:M) plati=D(i)>sum(abs(A(:,i)))-D(i); end //overeni pro sloupce ODD=isequal(plati, ones) //vysledek: 1 – ODD po sloupcich, 0 neni //normy norm(A) //vysledek: euklidovska norma norm(A,inf) //vysledek: radkova norma norm(A,1) //vysledek: sloupcova norma //kořeny polynomu roots([aN … a2 a1 a0]) // určení kořenů aN xN+…a2x2+a1x+a0 = 0

19

Prostá iterační metoda • Výpočet: x(k+1) = Ux(k)+v (pro U,v musí platit v=Qb, U+QA=E) • Podmínky konvergence: – Postačující : existuje alespoň 1 norma matice U, která je menší než 1. • Je-li splněna, lze odhadnout velikost chyby v kroku i:

|| U || || U ||i (i ) ( i 1) * (i ) || x  x || || x  x || nebo || x  x || || x (1)  x ( 0) || 1 || U || 1 || U || *

(i )

• Případně lze určit počet iterací k, abychom získali řešení s chybou menší než : || U || k   (1 || U ||) * (k ) || x  x || || x (1)  x ( 0) ||   k ln || U || ln (1) ( 0) 1 || U ||

|| x  x

||

– Nutná a postačující : spektrální poloměr matice U musí být menší než 1 (spektrální poloměr je největší z absolutních hodnot vlastních čísel matice). 20

Prostá iterační metoda-příklad  0.4  0.1 0  1     Dána soustava rovnic ve tvaru x=Ux+v, kde U   0.1 0.2 0.4 , v   2   0  0.4 0.4   3     a) Ověřte, že soustavu lze řešit prostou iterační metodou. Existuje norma matice U, která je menší než 1 (např. řádková: ||U||=0.8), je splněna postačující podmínka konvergence. Volte x(0)=(0,0,0)T a vypočtěte x(1) a x(2) prostou iterační metodou.

b)

x (1)

c)

U U x( 0) v x (1) v            0.4  0.1 0  0   1   1   0.4  0.1 0  1  1   1 .2         ( 2)          0.1 0.2 0.4  0    2    2 , x   0.1 0.2 0.4  2   2    3.7  0  0.4 0.4  0   3   3   0  0.4 0.4  3  3   3 .4              

Odhadněte velikost chyby x(2).

|| U ||  0.8 || x*  x ( 2) ||   || x ( 2)  x (1) ||  || x*  x ( 2) ||   1 || U ||  1  0.8

 0.2    0.8 * ( 2) 1.7  6.4  1.7  || x  x ||   0 . 2  0.4   

d) Určete počet iterací potřebných k výpočtu x(k) s chybou < 0.1. || x  x *

(k )

|| U || k 0.8k 0.1 0.2 0.02 (1) ( 0) ||   0.1  || x  x ||   0.1  .3  0.1  0.8k   k ln 0.8  ln  k  22.45 1 || U ||  1  0.8 3 3

Je potřeba 23 iterací.

21

Kontrola v MATLABu - prostá iterační metoda • U=[0.4 -0.1 0; 0.1 0.2 0.4; 0 -0.4 0.4] //zadani matice U (3x3) po radcich • max(abs(eig(A)))<1 // overeni konvergence • v=[1;2;3] // zadani vektoru – sloupce • xk=[0;0;0] //zadani x0 ////////////////// zpusob iterace po iteraci: • x_k1=U*xk+v //zobrazi se vysledek dalsi iterace • xk=x_k1; // strednik zpusobi nezobrazeni vysledku • x_k1=U*xk+v //pocitame dalsi iteraci • xk=x_k1; // … a tak dal… ////////////////// nebo pocitame urcity pocet iteraci: • pocet_iteraci=5 //kolik iteraci chceme vypocitat • for(i=1:1:pocet_iteraci) x_k1=U*xk+v; xk=x_k1;end //v xk bude posledni vypoctena ////// nebo pocitame dokud norma rozdilu 2 iteraci po sobe bude mensi nez cislo • cislo = 0.001 • x_k1=U*xk+v • while(norm(xk-x_k1)>cislo) xk=x_k1; xk_1=U*xk+v; end 22

Jacobiova metoda a a  • Iterační matice UJ     0 a a (formálně UJ =-D-1(L+U), vJ=D-1 b kde D,U,L jsou matice UJ A=L+D+U viz Gaussova – Seidelova metoda)

  a21  a  22   a31   a33      an1  a  nn

12

13

11

11

0 a32 a33  a  n2 ann





a23  a22 0

 



an 3 ann



a1n   b1     a11  a  11  a2 n   b2   a  a22    22  a3n , v J   b3    a33  a33        b    n  0   ann   

(aij jsou koeficienty původní matice soustavy A)

• Podmínky konvergence: – Postačující : • Matice A je ostře diagonálně dominantní • existuje norma matice UJ, která je menší než 1

– Nutná a postačující : • spektrální poloměr matice UJ musí být menší než 1 23

Jacobiova metoda - praktický výpočet • Z první rovnice vyjádříme 1. neznámou • Z druhé rovnice vyjádříme 2. neznámou • Ze třetí rovnice vyjádříme 3. neznámou …

x1( k 1)



x2( k 1)



x3( k 1)







xn( k 1)















1 b1  a12 x2k  a13 x3k    a1n xnk a11 1 b2  a21x1k  a23 x3k    a1n xnk a22 1 b3  a31x1k  a32 x2k    a1n xnk a33  1 bn  a11x1k  a12 x2k    a1n 1 xnk1 ann



• Postupně (v libovolném pořadí) vypočítáme nové, přesnější, hodnoty x(k+1). Do výrazů na pravé straně dosazujeme za x(k) hodnoty vypočítané na předcházejícím kroku. • Určení spektrálního poloměru iterační matice UJ z matice soustavy A: Vlastní čísla UJ matice určíme z rovnice: a13  a1n   a11 a12    a21 a22 a23  a2 n  det  a31 a32 a33  a3n   0        a  a a   a n 1 n 2 n 3 nn   24



Jacobiova metoda - příklad  4 1 0 1     Dána soustava rovnic Ax=b, kde A   2 4 1 , b    1  0 2 4 6     a) Ověřte, že danou soustavu lze řešit Jacobiovou metodou. b) Určete spektrální poloměr Jacobiovy matice. c) Volte x(0)=(0,0,0)T a určete x(1) a x(2) Jacobiovou metodou.

a) Matice A je ostře diagonálně dominantní, je splněna postačující podmínka 0  4 1 konvegrence Jacobiovy metody.   b) Vlastní čísla matice UJ určíme z výrazu det  2 4 1   0  0  2 4   

643  8  8  0  16 (42  1)  0  1  0, 2 

1 1 , 3   2 2

tedy (UJ)=0.5 25

Jacobiova metoda - příklad c) Z první rovnice vyjádříme x1, z druhé x2, ze třetí x3.

1 1  x2  4 Při výpočtu první iterace x(1)=(x1(1), x2(1), x3(1)) (levá strana) 1 x2   1  2 x1  x3  dosadíme do výrazů na pravé straně hodnoty x(0)=(0,0,0). 4 1 x3  6  2 x2  Dostáváme x(1)=(1/4,-1/4,3/2)T. 4 x1 

Při výpočtu druhé iterace x(2)=(x1(2), x2(2), x3(2)) (levá strana) dosadíme do výrazů na pravé straně hodnoty x(1). Dostáváme 1  1  5 x1( 2 )  1    4 4  16 1 1 3 3 x (22 )    1  2      4 4 2 4 1 1  13 x (32 )   6  2    4 4  8

26

Jacobiova metoda – kontrola v MATLABu A=[4 1 0; 2 4 1; 0 2 4]; b=[1;-1;6]; x0=[0;0;0]; ////vytvorime matici UJ a vektor vJ U=triu(A,1); L=tril(A,-1); D=A-U-L; UJ=inv(D)*(U+L)*(-1); vJ=inv(D)*b; //// vypocet iteraci jako prostou iteracni metodou x1=UJ*x0+vJ X2=UJ*x1+vJ ///... a tak dale rho=max(abs(eig(UJ))) ////spektralni polomer UJ vypocitany z matice A roots([64 0 -16 0 ]) rho=max(abs(roots ([64 0 -16 0 ]) ))

//spektralni polomer UJ

//urceni korenu 64*x3-16x=0

27

•

Gaussova-Seidelova metoda

Iterační matice UG=-(D+L)-1P (dokážeme se bez ní obejít!), vG=(D+L)-1b, kde

 a11 a12   a21 a22 A   a31 a32    a  n1 an 2  0 a12  0 0 U  0 0   0 0 

•

a13  a1n   b1     a23  a2 n   b2  a33  a3n , b   b3 ,        b  an 3  ann   n

a13  a1n   a11 0   a23  a2 n   0 a22 0 0  a3n , D   0         0 0 0  0  

0  0  0  0   0  0   a21 0 a33  0 , L   a31 a32        a 0  ann   n1 an 2

0 0 0 an 3

    

0  0 0   0 

Podmínky konvergence: – Postačující : • Matice A je ostře diagonálně dominantní • Matice A je symetrická a zároveň pozitivně definitní • existuje norma matice UG, která je menší než 1 – Nutná a postačující : • spektrální poloměr matice UG musí být menší než 1 28

•

Gaussova-Seidelova metoda - praktický výpočet 1  x  b  a x  a x  a x  Z první rovnice vyjádříme 1. neznámou a ( k 1) 1

11

•

( k 1) Z druhé rovnice vyjádříme 2. neznámou x2

•

Ze třetí rovnice vyjádříme 3. neznámou x3( k 1)

… •

1

 xn( k 1)

k 12 2

k 13 3



k 1n n



1 b2  a21x1k 1  a23 x3k    a1n xnk a22 1  b3  a31x1k 1  a32 x2k 1    a1n xnk a33   1  bn  a11x1k 1  a12 x2k 1    a1n 1 xnk11 ann 







Postupně (v pořadí od x1 k xn) vypočítáme nové, přesnější, hodnoty x(k+1). Do výrazů na pravé straně dosazujeme za xi hodnoty vypočtené na tomto kroku (pokud jsou k dispozici) nebo na předcházejícím kroku, tj.

při výpočtu xi(k+1) dosadíme x1(k+1) x2(k+1) …xi-1(k+1) a xi+1(k) … xn(k) . •



Určení spektrálního poloměru iterační matice UG z matice soustavy A: Vlastní čísla UG matice určíme z rovnice: a13   a11 a12

  a21 a22 det  a31 a32     a  n1 an 2

a23  a33   an 3 

a1n   a2 n  a3n   0    ann 

29

Gaussova-Seidelova metoda – příklad

 4 1 0 1     A   2 4 1 , b    1  0 2 4 6     a) Ověřte, že danou soustavu lze řešit Gaussovou-Seidelovou metodou. b) Určete spektrální poloměr Gaussovy-Seidelovy matice. c) Volte x(0)=(0,0,0)T a určete x(1) a x(2) Gaussovou-Seidelovou metodou. a) Matice A je ostře diagonálně dominantní, je splněna postačující podmínka konvegrence Gaussovy-Seidelovy metody. 0  4 1   b) Vlastní čísla matice UG určíme z výrazu det  2 4 1   0  0 2 4    Dána soustava rovnic Ax=b, kde

1 64  8  8  0  16 (4  1)  0  1  0, 2  0, 3  4 3

2

2

2

tedy (UG)=0.25 30

Gaussova-Seidelova metoda – příklad c) Z první rovnice vyjádříme x1, z druhé x2, ze třetí x3.

x1 

1 1  x2 , x2  1  1  2 x1  x3 , x3  1 6  2 x2  4 4 4

Výpočet první iterace:

( x2  83 )Výpočet

( x2  0 )

   1 x1(1)  1 4 ( x1  14 , x3  0 )    1 1  x (21)    1  2   0  4 4   3 ( x2  8 )  1 3 x (31)   6  2   4 8 

Změny vektoru x(0) na x(1) x ( 0)

 



druhé iterace:

 1  3 1 x1( 2 )  1   4 8   4 ( x1  13 12 , x3  12 67 )  3 1 11 27   x (22 )    1  2   8 4 32 16   2 7 ( x2  3 2 )  1  27  27 x3( 2 )   6  2   4 32  16



11 32







123 64

27 32

Změny vektoru x(1) na x(2)

 0  1 4   1 4   1 4   1 4   11 32   11 32   11 32                    0    0     3 8     3 8   x (1) x (1)    3 8     3 8     27 32     27 32   x ( 2)  27 16   27 16   27 16   123 64   0   0   0   27 16                  31

Gaussova-Seidelova metoda – kontrola v MATLABu A=[4 1 0; 2 4 1; 0 2 4]; b=[1;-1;6]; x0=[0;0;0]; ////vytvorime matici UGa vektor vG U=triu(A,1); L=tril(A,-1); D=A-U-L; UG=inv(D+L)*U*(-1); vG=inv(D+L)*b; //// vypocet iteraci jako prostou iteracni metodou x1=UG*x0+vG X2=UG*x1+vG ///... a tak dale rho=max(abs(eig(UJ))) ////spektralni polomer UG vypocitany z matice A roots([64 -16 0 0 ]) rho=max(abs(roots ([64 -16 0 0 ]) ))

//spektralni polomer UJ

//urceni korenu 64*x3-16x2=0

32

Řešení soustav lineárních rovnic gradientní metody (Podle F.Bubeník, M.Pultar, I.Pultarová – Matematické vzorce a metody, Vydavatelství ČVUT, 2001) • Metodu největšího spádu a metodu sdružených gradientů lze použít pro řešení soustav se symetrickou a pozitivně definitní maticí A. Řešení takové soustavy je jediným minimem kvadratické formy Q(x)=1/2xTAx-xTb. Volíme počáteční aproximaci x(0) . Aproximaci x(k+1) počítáme z aproximace x(k) přičtením (k) násobku vektoru p(k). Směr p(k) zvolíme a číslo (k) vypočítáme tak, aby aproximace x(k+1) minimalizovala kvadratickou formu Q na přímce x= x(k) + p(k). Metody se liší podle volby směru p(k). • Reziduum k-té aproximace je definováno vztahem r(k)=A x(k) – b • x (k+1)= x(k) + p(k) 33

Metoda největšího spádu • Směr p(k) se volí jako směr gradientu funkce Q.

x

( k 1)

x

( k )T

(k )

(k )

r r (k )  ( k )T ( k ) r r Ar

• Postačující podmínka konvergence : A je symetrická a pozitivně definitní. Výpočet 1. a 2. aproximace podrobněji: volíme x ( 0) vypočteme:

máme x (1) vypočteme:

reziduum r ( 0 )  Ax ( 0 )  b

reziduum r (1)  Ax (1)  b

skal. součin r ( 0 )T r ( 0) a součin r ( 0 )T  A  r ( 0) skal. součin r (1)T r (1) a součin r (1)T  A  r (1) (1)T (1) r ( 0 )T r ( 0 ) r r číslo   ( 0)T ( 0 ) číslo   (1)T (1) r Ar r Ar x (1)  x ( 0)  r ( 0 ) ( 2) (1) (1)

x

 x  r

34

Metoda největšího spádu - příklad Dána soustava rovnic Ax=b, kde

2 4 1  1     A   4 9 2 , b    3  1 2 5  4     

a) Ověřte, že danou soustavu lze řešit metodou největšího spádu. b) Volte x(0)=(0,0,0)T a určete x(1) metodou největšího spádu. a) Matice A je symetrická. Ověříme, zda je pozitivně definitní: a11 = 2 > 0, 18-16=2>0, det A = 90+8+8-9-8-80 = 9 > 0. Matice A je symetrická a pozitivně definitní, lze použít metodu největšího spádu. b) Výpočet x(1) : r(0)=Ax(0)-b= (1, 3, -4)T, r(0)Tr(0)=1  1+3  3+(-4)(-4) = 26, r(0)TAr(0) = (1, 3, -4)(2+12-4, 4+27-8, 1+6-20)T = (1, 3, -4)(10, 23, -13)T = 10+69+52=131, = 26/131 , x(1) = x(0) -  r(0) = -26/131 (1, 3, -4)T = (-0.198, -0.595, 0.794) T 35

Metoda sdružených gradientů • Směry p(k) se získají A-ortogonalizací posloupnosti vektorů reziduí r(k) . • Postačující podmínka konvergence : A je symetrická a pozitivně definitní. ( k 1)

 (k )

x  p Výpočet aproximací x : první aproximace – x(1) - analogicky metodě největšího spádu (k )

(k )

volíme x ( 0) ,

vypočteme: reziduum r ( 0)  Ax ( 0)  b, směr p ( 0)   r ( 0) ( 0 )T ( 0 ) r r číslo  ( 0)  ( 0)T ( 0) , aproximaci x (1)  x ( 0)   ( 0) p ( 0) p Ap druhá a další aproximace:

reziduum r (1)  r ( 0 )   ( 0 ) Ap ( 0 ) číslo 

(0)

směr p (1) číslo 

(1)

r (1)T r (1)  ( 0 )T ( 0 ) , r r  r (1)   ( 0 ) p ( 0 ) , r (1)T r (1)  (1)T (1) , p Ap

x ( 2 )  x (1)   (1) p (1)

reziduum r ( k 1)  r ( k )   ( k ) Ap ( k ) r ( k 1)T r ( k 1) číslo   ( k )T ( k ) , r r směr p ( k 1)  r ( k 1)   ( k ) p ( k ) , (k )

…

číslo 

( k 1)

r ( k 1)T r ( k 1)  ( k 1)T ( k 1) , p Ap

x ( k  2 )  x ( k 1)   ( k 1) p ( k 1) 36