Kapitola 3. SLAR - přímé metody

Kapitola 3. SLAR - pˇr´ım´ e metody

Formulace:

Je d´ana ˇctvercov´a matice A= [aij]ⁿ_i;j=1 a vektor prav´e strany b= [b1; b2; : : : ; bn]^T. Hled´ame vektor x = [x1; x2; : : : ; xn]^T tak, aby platilo

Ax= b;

rozeps´ano po sloˇzk´ach

a11x1+ a12x2+


+ a1nxn = b1

a21x1+ a22x2+


+ a2nxn = b2

...

an1x1+ an2x2+


+ annxn = bn

Pˇredpokl´ad´ame, ˇze je matice A regul´arn´ı (tj. soustava m´a pr´avˇe jedno ˇreˇsen´ı).

M´ame dva z´akladn´ı typy soustav:

 soustavy s obecnou matic´ı

 soustavy se speci´aln´ı matic´ı (symetrick´a, pozitivnˇe definitn´ı, ˇr´ıdk´a, p´asov´a apod.)

Pro prvn´ı skupinu se vˇetˇsinou pouˇz´ıvaj´ı pˇr´ım´e metody, pro druhou skupinu metody iteraˇcn´ı nebo speci´aln´ı modifikace pˇr´ım´ych metod.

Cramerovo pravidlo

nezn´am´a sloˇzka ˇreˇsen´ı xi= det A_i det A poˇcet operac´ı:

Je nutn´e vypoˇc´ıtat(n + 1) determinant˚u.

Pro v´ypoˇcet determinantu je tˇreba n! sˇc´ıt´an´ı a v kaˇzd´em sˇc´ıtanci je (n 1) n´asoben´ı.

Dost´av´ame:

(n + 1) [(n 1)n! + n!] = n(n + 1)!

napˇr: pron = 30, 10⁶ operac´ı za sekundu ! v´ypoˇcet trv´a 7; 82
10²¹ let

Idea dalˇs´ıch pˇr´ım´ych metod vych´azej´ı z faktu, ˇze soustavy

Ax= b a TAx = Tb , kde T je regul´arn´ı matice, maj´ı tot´eˇz ˇreˇsen´ı, tj. jsou ekvivalentn´ı.

Josef Danˇek 9.3.2o17 Touto transformac´ı lze z´ıskat troj´uheln´ıkovou soustavu

Ux= y : U= TA; y= Tb

pˇr: 2

64

u11 u12 u13

0 u22 u23

0 0 u33

3 75

2 64

x1

x2

x3

3 75=

2 64

y1

y2

y3

3 75

Troj´uheln´ıkovou soustavu lze velmi snadno ˇreˇsit zpˇetnou substituc´ı. Realizovan´y proces se naz´yv´a zpˇetn´y chod.

Gaussova eliminaˇcn´ı metoda

Ax= b rozeps´ano po sloˇzk´ach

2 64

a11 a12 a13 b1

a21 a22 a23 b2

a31 a32 a33 b3

3 75

Definujeme multiplik´atory m21 = a21

a₁₁ , m31 = a31

a₁₁ ˇr⁽¹⁾₁ = ˇr1 b⁽¹⁾₁ = b1

ˇr⁽¹⁾₂ = ˇr2+ m21ˇr₁ b⁽¹⁾₂ = b2+ m21b1

ˇr⁽¹⁾₃ = ˇr3+ m31ˇr₁ b⁽¹⁾₃ = b3+ m31b1

Z´ısk´ame novou soustavu : : : 1. f´aze eliminace

A⁽¹⁾x= b⁽¹⁾

2 66 66 4

a11 a12 a13 b1

0 a⁽¹⁾₂₂ a⁽¹⁾₂₃ b⁽¹⁾₂ 0 a⁽¹⁾₃₂ a⁽¹⁾₃₃ b⁽¹⁾₃

3 77 77 5

Definujeme multiplik´ator m₃₂ = a⁽¹⁾₃₂ a⁽¹⁾₂₂

ˇr⁽²⁾₁ = ˇr⁽¹⁾₁ b⁽²⁾₁ = b⁽¹⁾₁ ˇr⁽²⁾₂ = ˇr⁽¹⁾₂ b⁽²⁾₂ = b⁽¹⁾₂

ˇr⁽²⁾₃ = ˇr⁽¹⁾₃ + m₃₂ˇr⁽¹⁾₂ b⁽²⁾₃ = b⁽¹⁾₃ + m₃₂b⁽¹⁾₂

Z´ısk´ame novou soustavu : : : 2. f´aze eliminace

A⁽²⁾x= b⁽²⁾

2 66 66 4

a₁₁ a₁₂ a₁₃ b₁ 0 a⁽¹⁾₂₂ a⁽¹⁾₂₃ b⁽¹⁾₂ 0 0 a⁽²⁾₃₃ b⁽²⁾₃

3 77 77 5

Cel´y tento postup naz´yv´ame pˇr´ım´y chod. Troj´uheln´ıkovou soustavu ˇreˇs´ıme zpˇetn´ym chodem.

Efektivnost algoritmu GEM

Bereme v ´uvahu pouze operace n´asoben´ı a dˇelen´ı (poˇcet operac´ı sˇc´ıt´an´ı je pˇribliˇznˇe stejn´y).

(N . . . je ˇr´ad matice A)

• Celkem je N 1 f´az´ı eliminace. V K-t´e f´azi poˇc´ıt´ame N K multiplik´ator˚u (tj. N K dˇelen´ı)

N 1X

K=1(N K) = (N 1)N ^{N 1X}

K=1K = (N 1)N 1

2(N 1)N = 1

2(N 1)N

• Kaˇzd´ym multiplik´atorem vyn´asob´ıme(N K + 1) prvk˚u rozˇs´ıˇren´e matice (jeden rozˇs´ıˇren´y ˇr´adek), tj. (N K)(N K + 1) v K-t´e f´azi

N 1X

K=1

(N K)(N K + 1) =^{N 1X}

K=1

h(N² + N) K(2N + 1) + K²ⁱ=

= (N 1)(N²+ N) 1

2N(N 1)(2N + 1) +1

6(N 1)N(2N 1) =

= N³ N²+ N² N N³+1

2N²+ 1 2N +1

3N³ 1

2N² +1 6N =

= 1

3N³ 1 3N

• Zpˇetn´y chod

1 + ( 1_|{z}

dˇelen´ı

+ 1_|{z}

n´asoben´ı

) + (1 + 2) +


+ (1 + N 1) = ^XN

K=1K = 1

2N(N + 1)

Celkem

1

2N(N 1)

| {z }

v´ypoˇcet multiplik´ator˚u

+1

3N³ 1 3N

| {z }

pˇr´ım´y chod

+1

2N(N + 1)

| {z }

zpˇetn´y chod

= 1

3N³+ N² 1 3N

Pˇr´ıklad 1

Reˇste soustavu rovnicˇ

x + 2y + 3z = 14 2x + 4y + 5z = 25 7x + 8y + 9z = 50

; tj.

2 64

1 2 3 2 4 5 7 8 9

3 75

2 64

x y z

3 75=

2 64

14 25 50

3 75

Reˇˇ sen´ı

Josef Danˇek 9.3.2o17 Pro z´apis budeme pouˇz´ıvat tvar matice rozsrene:

2 66 4

1 2 3 14 2 4 5 25 7 8 9 50

3 77 5

=
( ²₁) ^

+ =
( ₁⁷) ¹

A+

2 66 64

1 2 3 14

0 0 1 3

0 6 12 48

3 77

75 =
( ⁶₀) ^

+

!!!

 Algoritmus Gaussovy eliminaˇcn´ı metody pro tento pˇr´ıklad nen´ı realizovateln´y.

 Snadno se pˇresvˇedˇc´ıme, ˇze m´a dan´a soustava ˇreˇsen´ı

x = 1; y = 2 a z = 3:

Gaussova eliminaˇcn´ı metoda ale selhala.

Ot´azka 1: Pro jak´e matice A m´a soustava Ax = b pr´avˇe jedno ˇreˇsen´ı?

! matice A mus´ı b´yt regul´arn´ı, tj.

vˇsechna vlastn´ı ˇc´ısla mus´ı b´yt r˚uzn´a od nuly

jinak ˇreˇceno ˇr´adky matice A mus´ı b´yt line´arnˇe nez´avisl´e jinak ˇreˇceno sloupce matice A mus´ı b´yt line´arnˇe nez´avisl´e jinak ˇreˇceno det A 6= 0

Pozn´amka: Vlastn´ı ˇc´ıslo matice A je ˇc´ıslo  splˇnuj´ıc´ı rovnici Av = v, kde v je vlastn´ı vektor odpov´ıdaj´ıc´ı vlastn´ımu ˇc´ıslu . ˇC´ıslo  tedy urˇcit´ym zp˚usobem charakterizuje matici A.

Ot´azka 2: Pro jak´e matice A je algoritmus Gaussovy eliminaˇcn´ı metody realizovateln´y?

! Vˇeta: Je-li matice A ostˇre diagon´alnˇe dominantn´ı, pak je algoritmus GEM realizovateln´y.

Pozn´amka: Matice A= [aij]i;j=1;:::;n je ostˇre diagon´alnˇe dominantn´ı, plat´ı-li ja_iij > ^Xn

j=1;j6=i

ja_ijj;

tj. absolutn´ı hodnota diagon´aln´ıho prvku je vˇetˇs´ı neˇz souˇcet absolutn´ıch hodnot ostatn´ıch prvk˚u v ˇr´adku.

Napˇr.:

A=

2 64

7 1 2 1 4 1 2 4 9

3 75:

! Vˇeta: Je-li matice A symetrick´a a pozitivnˇe definitn´ı, pak je algoritmus GEM realizovateln´y.

Pozn´amka: Matice A je symetrick´a, plat´ı-li pro jej´ı prvky

aij = aji 8i; j = 1; 2; : : : ; n:

Pozn´amka: Matice A je pozitivnˇe definitn´ı, m´a-li vˇsechna vlastn´ı ˇc´ısla kladn´a

nebo jinak ˇreˇceno 8x 6= o : x^TAx> 0.

Pozn´amka: Pro soustavu s matic´ı, kter´a splˇnuje pˇredpoklady nˇekter´e z uveden´ych vˇet, je moˇzn´e dopˇredu ˇr´ıci, ˇze p˚ujde ˇreˇsit pomoc´ı Gaussovy eliminaˇcn´ı metody. Obr´acenˇe to ovˇsem neplat´ı, tj. nen´ı-li napˇr.

matice soustavy ostˇre diagon´alnˇe dominantn´ı, jeˇstˇe to obecnˇe neznamen´a, ˇze nep˚ujde pomoc´ı Gaussovy eliminaˇcn´ı metody ˇreˇsit.

Pozn´amka: Abychom zaruˇcili, ˇze soustava p˚ujde vyˇreˇsit pro libovolnou regul´arn´ı matici, mus´ıme algo- ritmus Gaussovy eliminaˇcn´ı metody upravit. Zavedeme tzv. v´ybˇer hlavn´ıho prvku (pivotaci).

Pozn´amka: Pivot (hlavn´ı prvek) . . . prvn´ı nenulov´y prvek v dan´em ˇr´adku matice.

Pˇr´ıklad 2

Pomoc´ıGEM se sloupcovou pivotac´ıvyˇreˇste soustavu rovnic z Pˇr´ıkladu 1, kde selhala klasick´a GEM, tj. ˇreˇs´ıme soustavu Ax= b, kde

A=

2 64

1 2 3 2 4 5 7 8 9

3

75 a b =

2 64

14 25 50

3 75:

Reˇˇ sen´ı

1. sloupec

#

vymˇeˇn

2 4

!

!

0 BB B@

1 2 3 14 2 4 5 25 7 8 9 50

1 CC CA ;

0 BB B@

7 8 9 50 2 4 5 25 1 2 3 14

1 CC CA

=
( ²₇) ^

+ =
( ¹₇) ¹

A+

2. sloupec

#

Josef Danˇek 9.3.2o17

nen´ı tˇreba mˇenit

!

0 BB BB

@

7 8 9 50 0 ¹²₇ ¹⁷₇ ⁷⁵₇ 0 ⁶₇ ¹²₇ ⁴⁸₇

1 CC CC A

=
( 12⁶⁷

7 ) = ¹₂ ^! +

0 BB BB

@

7 8 9 50 0 ¹²₇ ¹⁷₇ ⁷⁵₇ 0 0 ¹₂ ³₂

1 CC CC A

=) x=

2 64

1 2 3

3 75

Pozn´amky:

• Pˇri sloupcov´e pivotaci jsme postupnˇe v kaˇzd´em sloupci (resp. jeho ˇc´asti pod diagon´alou vˇcetnˇe) vyb´ırali ˇc´ıslo, kter´e bylo maxim´aln´ı v absolutn´ı hodnotˇe a v pˇr´ıpadˇe, ˇze toto ˇc´ıslo neleˇzelo na diagon´ale, vymˇenili jsme pˇr´ısluˇsn´e 2 rovnice. D´ale jsme pokraˇcovali jako v GEM bez pivotace, tj. nulovali jsme koeficienty pod diagon´alou.

• Sloupcov´a pivotace nen´ı jedin´a moˇznost. Podobnˇe m˚uˇzeme vyb´ırat i maxim´aln´ı prvek v absolutn´ı hodnotˇe z pˇr´ısluˇsn´eho ˇr´adku (resp. jeho ˇc´asti) a pot´e vymˇenit pˇr´ısluˇsn´e sloupce. Pozor! Je ovˇsem tˇreba zamˇenit i pˇr´ısluˇsn´e sloˇzky ˇreˇsen´ı x. V tomto pˇr´ıpadˇe hovoˇr´ıme oˇr´adkov´e pivotaci.

• Dalˇs´ı moˇznost´ı je vyb´ırat maxim´aln´ı prvek v absolutn´ı hodnotˇe z cel´e matice A (resp. pˇr´ısluˇsn´e podma- tice). V tomto pˇr´ıpadˇe hovoˇr´ıme o ´upln´e pivotaci. Opˇet je tˇreba m´ıt na pamˇeti, ˇze je tˇreba zamˇenit sloˇzky ve vektoru ˇreˇsen´ı. Nev´yhodou ´upln´e pivotace je pomalejˇs´ı v´ypoˇcet nebot’ hlavn´ı prvek vyhled´av´ame z cel´e dosud neupraven´e ˇc´asti.

• Libovolnou pivotac´ı dos´ahneme realizovatelnosti GEM pro libovolnou regul´arn´ı matici A.

Metoda LU-rozkladu

Opˇet uvaˇzujeme regul´arn´ı matici A ˇr´adu N. Matici A lze rozloˇzit na souˇcin A = LU , kde L je doln´ı troj´uheln´ıkov´a matice ˇr´adu N a U je horn´ı troj´uheln´ıkov´a matice ˇr´adu N.

Napˇr:

2 64

a₁₁ a₁₂ a₁₃ a21 a22 a23

a31 a32 a33

3 75=

2 64

l₁₁ 0 0 l21 l22 0 l31 l32 l33

3 75

2 64

u₁₁ u₁₂ u₁₃ 0 u22 u23

0 0 u33

3 75

Tento rozklad nen´ı d´an jednoznaˇcnˇe (12 nezn´am´ych a 9 podm´ınek), jednoznaˇcnosti dos´ahneme napˇr. t´ım, ˇze poloˇz´ıme lii = 1; i = 1; 2; : : : ; N .

Algoritmus: (viz skripta)

je odvozen z postupn´eho n´asoben´ı ˇr´adk˚u matice L a sloupc˚u matice U

(1; 1) u11= a11 (2; 1) a21 = l21u11 (3; 1) a31 = l31u11

(1; 2) u12= a12 (2; 2) a22 = l21u12+ u22 (3; 2) a32 = l31u12+ l32u22

(1; 3) u13= a13 (2; 3) a23 = l21u13+ u23 (3; 3) a33 = l31u13+ l32u23+ u33

Reˇsen´ı soustavy Axˇ = b metodou LU-rozkladu:

1. Realizace LU-rozkladu: A= LU

2. ˇReˇsen´ı troj´uheln´ıkov´e soustavy: Ly= b L Ux_|{z}

y

= b 3. ˇReˇsen´ı troj´uheln´ıkov´e soustavy: Ux= y

Souvislost GEM a metody LU-rozkladu

Gaussovu eliminaci lze popsat pomoc´ı n´asoben´ı regul´arn´ımi maticemi.

Prvn´ı f´azi pop´ıˇseme takto : : : A⁽¹⁾ = M1A , kde M₁ =

2 66 66 66 66 66 4

1 m21 1

m31 1

m₄₁ 1

... . ..

m_n1 1

3 77 77 77 77 77 5

2 66 66 66 64

a₁₁ a₁₂ a₁₃ : : : a_1n

m21a11+ a21 m21a12+ a22 m21a13+ a23 : : : m21a1n+ a2n

m31a11+ a31 m31a12+ a32 m31a13+ a33 : : : m31a1n+ a3n

... ... ... ...

mn1a11+ an1 mn1a12+ an2 mn1a13+ an3 : : : mn1a1n+ ann

3 77 77 77 75

=

=

2 66 66 66 66 66 4

1 m₂₁ 1

m₃₁ 1

m41 1

... . ..

mn1 1

3 77 77 77 77 77 5

2 66 66 66 64

a11 a12 a13 : : : a1n

a21 a22 a23 : : : a2n

a31 a32 a33 : : : a3n

... ... ... ... an1 an2 an3 : : : ann

3 77 77 77 75

Druhou f´azi pop´ıˇseme takto : : : A⁽²⁾ = M2A⁽¹⁾ , kde M₂ =

2 66 66 66 66 66 4

1 1 m32 1

m42 1

... . ..

mn2 1

3 77 77 77 77 77 5

...

Josef Danˇek 9.3.2o17

Nakonec(n 1) f´azi pop´ıˇseme : : : A^{(n 1)} = Mn 1A^{(n 2)} , kde M_{n 1} =

2 66 66 66 66 66 4

1 1

1 . ..

1 m_{n;n 1} 1

3 77 77 77 77 77 5

Dostali jsme horn´ı troj´uheln´ıkovou matici, oznaˇc´ıme ji napˇr. V V= A^{(n 1)} = M_{n 1}M_{n 2}M_{n 3}: : : M₂M₁

| {z }

ozn.M

A ) A = M ¹V

M ¹ = M₁¹M₂¹: : : M_{n 1}¹

Jak vypad´a napˇr. M₂¹ ?

M₂ =

2 66 66 66 66 66 4

1 1 m32 1

m₄₂ 1

... . ..

mn2 1

3 77 77 77 77 77 5

! M₂¹ =

2 66 66 66 66 66 66 66 66 4

1 1

m32 1

m42 1

... . ..

mn2 1

3 77 77 77 77 77 77 77 77 5

protoˇze po vyn´asoben´ı:

• i-t´y ˇr´adek  j-t´y sloupec (j 6= i) – bud’ m₄₂
1 + 1
( m₄₂) = 0

– nebo m42
0 + 0
1 + 1
0 = 0

• i-t´y ˇr´adek  i-t´y sloupec – m42
0 + 1
1 = 1

tj. M₂
M₂¹ = I

Jak vypad´a M ¹ ?

M ¹ =

2 66 66 66 64

1

m₂₁ 1

m₃₁ m₃₂ 1

... . ..

m_n1 m_n2 m_n3 : : : 1

3 77 77 77 75

Pˇr.:

2 64

1 0 0

m21 1 0 m31 0 1

3 75

2 64

1 0 0

0 1 0

0 m32 1

3 75=

2 64

1 0 0

m21 1 0

m31 m32 1

3 75

Plat´ı

A= M_{| {z }}¹

()

V_|{z}

()

(*) doln´ı troj´uheln´ıkov´a s 1 na diagon´ale (**) horn´ı troj´uheln´ıkov´a

) jedn´a se o LU-rozklad (rozklad je jednoznaˇcn´y) L= M ¹ a U= V

V´ypoˇcet determinant˚u

1. Uˇzit´ı GEM

U= M_{N 1}M_{N 2}: : : M₂M₁A

det U = det M| {zN 1}

=1

det MN 2

| {z }

=1

: : : det M| {z }2

=1

det M1

| {z }

=1

det A

det A = det U =^YN

i=1uii

2. Uˇzit´ı LU-rozkladu

det A = det L det U = ^YN

i=1

u_ii

V´ypoˇcet inverzn´ı matice 1. Uˇzit´ı GEM

A X= I (maticov´a soustava) X= A ¹

Josef Danˇek 9.3.2o17 2. Uˇzit´ı LU-rozkladu

A= L U A ¹ = U ¹L ¹

Numerick´e aspekty GEM a metody LU-rozkladu

Pˇri numerick´e realizaci nevypoˇcteme pˇresnˇe matice L a U, ale pˇribliˇzn´e matice L a^e U.^f Teoreticky plat´ı A= LU.

Oznaˇc´ıme ^fA=L^efU . . . dopoˇcteno pro z´ıskan´e matice L a^{e f}U.

Budeme zkoumat rozd´ıl ^fA A.

Oznaˇcme E a F matice chyb takov´e, ˇze plat´ı:

Le = L + E; ^fU= U + F

Potom:

fA A=L^efU LU= (L + E)(U + F) LU = EU + LF + EF

Odtud plyne z´avˇer: Pokud jsou multiplik´atory v absolutn´ı hodnotˇe velk´e, pak prvky L jsou v absolutn´ı hodnotˇe velk´e ) chyba m˚uˇze b´yt velk´a. Toto je jeden z d˚uvod˚u realizace pivotace.

Pˇr´ım´e metody pro soustavy se speci´aln´ı matic´ı

Uvaˇzujeme matice:

• symetrick´a

• symetrick´a a pozitivnˇe definitn´ı

• diagon´alnˇe dominantn´ı

• p´asov´a

Plat´ı: Je-li matice A symetrick´a a A^(k) jsou matice z´ıskan´e GEM v z´akladn´ı verzi, pak podmatice A^(k) jsou tak´e symetrick´e.

pˇr:

2 66 66 4

a b c b d e c e f

3 77 77 5

=
( _a^b) ^

+ =
( ^c_a) ¹

A+

2 66 66 66 4

a b c

0 d ^b_a² e ^bc_a 0 e ^bc_a f ^c_a²

3 77 77 77 5

Lze pak pouˇz´ıt symetrickou verzi GEM a LU-rozkladu.

Je-li matice A nav´ıc pozitivnˇe definitn´ı, pak lze realizovat Cholesk´eho metodu rozkladu A = U^TU

Pozn´amka: V algoritmu je potˇreba realizovat v´ypoˇcet odmocnin a to lze pouze pro pozitivnˇe definitn´ı matice.

pˇr: 2

64

a11 a12 a13

a₂₁ a₂₂ a₂₃ a₃₁ a₃₂ a₃₃

3 75=

2 64

u11 0 0 u₂₁ u₂₂ 0 u₃₁ u₃₂ u₃₃

3 75

2 64

u11 u12 u13

0 u₂₂ u₂₃ 0 0 u₃₃

3 75

(1; 1) : a₁₁= u²₁₁ (2; 1) : a₂₁ = u₁₂
u₁₁ (3; 1) : a₃₁ = u₁₃
u₁₁

(1; 2) : a₁₂= u₁₁
u₁₂ (2; 2) : a₂₂ = u²₁₂+ u²₂₂ (3; 2) : a₃₂ = u₁₃
u₁₂+ u₂₃
u₂₂ (1; 3) : a13= u11
u13 (2; 3) : a23 = u12
u13+ u22
u23 (3; 3) : a33 = u²₁₃+ u²₂₃+ u²₃₃

Metoda LU-rozkladu pro p´asov´e matice

Uvaˇzujeme matici A takovou, ˇze a_ij = 0 , kdyˇz ji jj > p (ˇs´ıˇrka p´asu je 2p + 1).

Pokud lze realizovat LU-rozklad, pak

lij = 0, kdyˇz j > i a j < i p, u_ij = 0, kdyˇz j < i a j > i + p.

Josef Danˇek 9.3.2o17 Pozn´amka: V pˇr´ıpadˇe obecn´e matice vˇsak nelze ˇcekat, ˇze matice L a U bude m´ıt nulov´y prvek v t´eˇze pozici jako jej mˇela matice A.

Pro p´asov´e, symetrick´e, pozitivnˇe definitn´ı matice pouˇz´ıv´ame speci´aln´ı verzi Cholesk´eho rozkladu.

Metoda faktorizace pro tˇr´ıdiagon´aln´ı matici

Uvaˇzujeme soustavun + 1 line´arn´ıch algebraick´ych rovnic A
Y = F ve tvaru:

0 BB BB BB BB BB BB

@

c0 b0 f0

a1 c1 b1 f1

a₂ c₂ b₂ f₂

a₃ c₃ b₃ f₃

. .. . .. . .. ... a_{n 1} c_n+1 b_{n 1} f_{n 1}

an cn fn

1 CC CC CC CC CC CC A

Oznaˇc´ıme

1 = b0

c₀ ; 1 = f0

c₀

Prvn´ı dvˇe rovnice (prvn´ı rovnice je vydˇelen´a c0) lze ps´at ve tvaru:

y_{0 1}y₁ = ₁

a1y0 + c1y1 b1y2 = f1

=
a₁ ^ +

(c₁ a_11) y₁ b₁y₂ = f₁ + a_11

Pˇrep´ıˇseme (rovnice vydˇelen´a koeficientem u y₁) na tvar:

y_{1 2}y₂ = ₂ kde

2 = b₁

c1 a11 ; 2 = f₁+ a_11 c1 a11

Po zobecnˇen´ı:

• PˇR´IM´Y CHOD

₁ = b0

c_{0 1} = f0

c₀

_i+1 = b_i

ci ai
i i = 1; 2; : : : ; n 1 _i+1 = f_i+a_i
i

ci ai
i i = 1; 2; : : : ; n

• ZPˇETN´Y CHOD

Pro posledn´ı dvˇe rovnice z´ısk´ame:

yn 1 nyn = n

anyn 1 + cnyn | = fn

=
a_n ^ +

(cn ann) yn = fn + ann

| ve druh´e rovnici jiˇz nen´ı ˇclen bnyn+1

Vyj´adˇr´ıme posledn´ı sloˇzku ˇreˇsen´ı:

y_n = f_n+ a_nn

cn ann =: _n+1

Zpˇetnˇe dosazujeme:

yi 1= i+i
yi i = n; n 1; : : : ; 1 ()

Efektivnost algoritmu

dˇelen´ı 2n + 1 (1 + n 1 + 1 + n)

n´asoben´ı
3n (n + n + n)

sˇc´ıt´an´ı a odˇc´ıt´an´ı  3n (n + n + n)

celkem (8n + 1) operac´ı

Pozn´amka: Pokud budeme ˇreˇsit tuto soustavu pro r˚uzn´e prav´e strany F, nemus´ıme jiˇz znovu vyjadˇrovat koeficienty _i, protoˇze nez´avis´ı na F. Staˇc´ı tedy pˇrepoˇc´ıtat _i.

Zat´ım jsme neuvedli pˇredpoklady pro metodu faktorizace

• je tˇreba zajistit, aby jmenovatel ci aii byl nenulov´y pro i = 1; 2; : : : ; n

• yi se urˇcuje z rekurentn´ı formule (), pˇritom m˚uˇze doj´ıt k akumulaci zaokrouhlovac´ıch chyb.

Necht’ i, i jsou dokonce pˇresnˇe vypoˇc´ıtan´e a necht’ m´ame yen = yn+ "n (s chybou"n).

Potom postupnˇe vypoˇc´ıt´ame podle ()

ye_{i 1} = _i+ _iy^e_i ; i = n; n 1; : : : ; 1

Josef Danˇek 9.3.2o17 Oznaˇc´ıme-li "_i=y^e_i y_i chybu, bude jistˇe splˇnovat homogenn´ı rovnici

"_{i 1} = _i"_i ; i = n; n 1; : : : ; 1 (protoˇze pˇresn´e hodnotyyi splˇnuj´ıyi 1 = i+ iyi)

) Pokud by byly koeficienty jij > 1, dojde k velk´emu n´ar˚ustu chyby "0 !!!

Pro jij  1; i = 1; 2; :::; n

| {z }

()

je algoritmus stabiln´ı.

Postaˇcuj´ıc´ı podm´ınky pro zajiˇstˇen´ı(): Matice soustavy je ostˇre diagon´alnˇe dominantn´ı.

D˚ukaz viz literatura.

Pˇr´ıklady aplikac´ı, kter´e vedou na soustavu s tˇr´ıdiagon´aln´ı matic´ı

1. ˇReˇsen´ı okrajov´e ´ulohy

(k(x)u⁰(x))⁰ q(x)u(x) = f(x) x 2 (0; l) u(0) = 0

u(l) = 0 k(x) > 0 q(x)  0

x na h0; li diskretizujeme : : : s´ıt’x_i

p˚uvodn´ı ´ulohu nahrad´ıme ´ulohou s diferenˇcn´ı rovnic´ı a pouˇzijeme vzorec pro pomˇernou diferenci

2. Diferenˇcn´ı sch´emata pro rovnici veden´ı tepla

@u

@t = @²u

@x² x 2 (0; l) ; t > 0 u(0; t) = ₁(t)

u(l; t) = 2(t) u(x; 0) = u0(x)

3. Soustava pro v´ypoˇcet koeficient˚u kubick´eho spline (aproximace funkce) . . . budeme prob´ırat

Podm´ınˇenost ´ulohy ˇreˇsit SLAR

Uvaˇzujeme opˇet soustavu Ax= b, A . . . n  n regul´arn´ı, b 2 Rⁿ. Oznaˇcen´ı:

A . . . mal´a zmˇena matice A

b . . . mal´a zmˇena vektoru b

x . . . odpov´ıdaj´ıc´ı zmˇena vektoru nezn´am´ych x^ . . . pˇresn´e ˇreˇsen´ı soustavy Ax= b

Plat´ı:

(A +
A)(x^+
x) = b +
b 1. Uvaˇzujme situaci
A = 0, tj. A je zad´ana pˇresnˇe

Ot´azka: Jakou zmˇenu ˇreˇsen´ı vyvol´a zmˇena prav´e strany?

Ax^+ A
x = b +
b A
x =
b

x = A ¹
b Z vlastnost´ı maticov´e normy plyne:

Josef Danˇek 9.3.2o17

Ax^ = b ) kbk  kAk
kx^k ) 1

kx^k  kAk kbk

x = A ¹
b ) k
xk  kA ¹k
k
bk

k
xk

kx^k  kA ¹k
k
bk
kAk kbk

C_p =

k
xkkx^k k
bkkbk

 kA ¹k
kAk

2. Pˇr´ıpad, kdy
b = 0, tj. b je zad´ana pˇresnˇe

(A +
A)(x^+
x) = b

Ax^+
Ax^+ A
x +
A
x = b

A
x =
A(x^+
x)

x = A ¹
A(x^+
x)

k
xk  kA ¹k
k
Ak
kx^+
xk
kAk kAk k
xk

kx^+
xk  kA_| ¹k
kAk_{{z }}

Cp

k
Ak kAk

C_p=

kx^k
xk+
xk k
AkkAk

 kA ¹k
kAk

3. Rozmyslete obecn´y pˇr´ıpad
b 6= 0,
A 6= 0 viz skripta (D.cv.)

Pozn´amka: Pro symetrick´e matice je ˇc´ıslo podm´ınˇenosti pod´ıl nejvˇetˇs´ı a nejmenˇs´ı absolutn´ı hodnoty vlastn´ıho ˇc´ısla.

C_p= kA_{| {z }}¹k

min1

kAk_{| {z }}

max

= j_maxj jminj

(Pro symetrick´e matice plat´ı, ˇze maj´ı re´aln´a vlastn´ı ˇc´ısla.)

Geometrick´a interpretace - dobˇre podm´ınˇen´a ´uloha (2D)

x + y = 3

x y = 1 ) y = 3 x

y = x 1

9=

; ˇreˇsen´ı x^ = 2 y^ = 1

A=

"

1 1

1 1

#

; b =

"

3 1

#

Dobˇre podm´ınˇen´a ´uloha - mal´a relativn´ı zmˇena vstupn´ıch dat vyvol´a malou relativn´ı zmˇenu v´ystupn´ıch dat.

Cp= kA ¹k
kAk

A ¹ =

"

0; 5 0; 5 0; 5 0; 5

#

= 0; 5 A

C_p= 1
2 =2 (ˇr´adkov´a, sloupcov´a norma); C_p= ^p1₂p

2 =1 (spektr´aln´ı norma)

Josef Danˇek 9.3.2o17 A^TA =

"

2 0 0 2

#

kAkSP =p

2; kA ¹kSP = ^p1₂

!

1. zmˇena prav´e strany (zmˇena matice
A = 0)

b =

"

0; 1 0; 2

#

; b +
b =

"

3; 1 1; 2

#

zmˇenˇen´a soustava y = 3; 1 x; y = x 1; 2

2. zmˇena matice soustavy (zmˇena prav´e strany
b = 0)

A =

"

0 0; 1 0; 2 0

#

; A +
A =

"

1 1; 1 1; 2 1

#

zmˇenˇen´a soustava y = _1;1¹ (3 x); y = 1; 2 x 1

Geometrick´a interpretace - ˇspatnˇe podm´ınˇen´a ´uloha (2D)

x + y = 2

x + 1; 1 y = 2; 1 ) y = 2 x y = _1;1¹ (2; 1 x)

9=

; ˇreˇsen´ı x^ = 1 y^ = 1

Josef Danˇek 9.3.2o17

A=

"

1 1

1 1; 1

#

; b =

"

2 2; 1

#

Spatnˇˇ e podm´ınˇen´a ´uloha - mal´a relativn´ı zmˇena vstupn´ıch dat vyvol´a velkou relativn´ı zmˇenu v´ystupn´ıch dat.

C_p= kA ¹k
kAk

A ¹ =

"

11 10

10 10

#

Cp= 2; 1
21 =44,1(ˇr´adkov´a, sloupcov´a norma); Cp= 2; 0512
20; 5125 = 42,07 (spektr´aln´ı norma)

1. zmˇena prav´e strany (zmˇena matice
A = 0)

b =

"

0; 1 0; 2

#

; b +
b =

"

2; 1 2; 3

#

zmˇenˇen´a soustava y = 2; 1 x; y = _1;1¹ (2; 3 x)

2. zmˇena matice soustavy (zmˇena prav´e strany
b = 0)

A =

"

0 0

0 0; 05

#

; A +
A =

"

1 1

1 1; 05

#

zmˇenˇen´a soustava y = 2 x;y = _1;05¹ (2; 1 x)

Pozn´amka: Jin´y v´ypoˇcet ˇc´ısla podm´ınˇenosti (prakticky)

Ax= b zmˇena na vstupu : : :
A

vyvol´a zmˇenu na v´ystupu : : :
x

(A +
A)(x +
x) = b

Cp=

k
xkkxk k
AkkAk

Pˇr´ıklad

A=

"

50 100 50 101

#

; b =

"

0 1

#

; x=

"

2 1

#

zmˇena matice soustavy

A =

"

0 0; 1 0; 2 0

#

) ^fA= A +
A =

"

50 99; 9 50; 2 101

#

e

x=^fA ¹b=

"

2; 8527 1; 4278

#

Josef Danˇek 9.3.2o17

x = x x^e =

"

2 1

# "

2; 8527 1; 4278

#

=

"

0; 8527 0; 4278

#

k:k 1.vektorov´a / sloupcov´a maximov´a / ˇr´adkov´a euklidovsk´a / spektr´aln´ı

x =

"

0; 8527 0; 4278

#

1,2805 0,8527 0,9539

x=

"

2 1

#

3 2 2,2361

A =

"

0 0; 1 0; 2 0

#

0,2 0,2 0,2

A=

"

50 100 50 101

#

201 151 158,7479

Cp 428,97 321,89 338,6

kxk₁ =^X

i

jx_ij ; kxk₁= max

i jx_ij ; kxk₂ =^sX

i

x²_i

kAkS = max

k

X

i

jaikj ; kAkR= max

i

X

k

jaikj ; kAkSP = max

k _k¹²(A^HA) (Pˇri pouˇzit´ı r˚uzn´ych norem dost´av´ame obecnˇe r˚uzn´a ˇc´ısla podm´ınˇenosti.)

Cp= kA ¹k
kAk

A ¹ =

"

2; 02 2

1 1

#

Cp= 151
4; 02 = 607; 02 (ˇr´adkov´a k:k);

Cp = 201
3; 02 = 607; 02 (sloupcov´a k:k);

Cp = 158; 7479
3; 1750 := 504; 03 (spektr´aln´ı k:k)

Pozn´amky k podm´ınˇenosti

• St´ale pˇredpokl´ad´ame, ˇze A je regul´arn´ı matice.

Soustava Ax= b m´a potom pr´avˇe jedno ˇreˇsen´ı.

Pˇredpokl´adejme, ˇze je matice A normalizov´ana tak, ˇze jej´ı max. prvek v absolutn´ı hodnotˇe je roven1.

Je-li soustava ˇspatnˇe podm´ınˇen´a, potom matice A ¹ mus´ı obsahovat velk´e prvky.

napˇr:

A=

"

1 1

1 0; 99999

#

) A ¹ =

"

99999 100000 100000 100000

#

To odpov´ıd´a faktu, ˇze ˇc´ıslo podm´ınˇenosti je velk´e C_p= kAk
kA ¹k (norma kA ¹k je velk´a).

• Rozbor chyb

Ax= b , x^ . . . pˇresn´e, x_c : : : vypoˇcten´e

) chybu m˚uˇzeme mˇeˇrit pomoc´ı rezidua r = Ax_c b (pro pˇresn´e ˇreˇsen´ı je r= Ax^ b= 0)

Plat´ı: Je-li x_c bl´ızko x^ ) reziduum je mal´e. Bohuˇzel to neplat´ı obr´acenˇe !!!

r = A(xc x^)

x_c x^ = A ¹r

Pro ˇspatnˇe podm´ınˇenou ´ulohu obsahuje A ¹ velk´e prvky, kter´e i pro mal´e hodnoty r mohou znamenat velkou chybu.

pˇr:

A=

"

1 1

1 0; 99999

#

; b =

"

2 1; 99999

#

pˇresn´e ˇreˇsen´ı: x^ = [1; 1]^T

vypoˇcten´e m˚uˇze b´yt: x_c = [ 98; 100]^T

r = Axc b=

"

2 1; 999

# "

2 1; 99999

#

=

"

0 0; 00099

#

) je lepˇs´ı pokouˇset se odhadovat kxc x^k.

Bohuˇzel se v odhadech vˇzdy vyskytuje norma kA ¹k !!! Podrobnˇeji viz literatura.