Choleskyho rozklad

Download Report

Transcript Choleskyho rozklad

Metóda Konečných Prvkov
vo výrobných technológiach
prednáška č. 2
Obsah prednášky
• Vektor, matica, tenzor
• Základné matematické operácie
• Spôsoby riešenia sústavy lineárnych algebraických rovníc
– Príklady
• Spôsoby riešenia sústavy nelineárnych algebraických
rovníc
• Riešiče používané programom ANSYS
prednáška č.2 - 2/59
Základná terminológia
Tenzor
•
•
je fyzikálna veličina nezávislá od aktuálne definovaného súradnicového
systému
je určená stupňom a usporiadaním
Matica
•
tenzor 2-ého rádu (stupňa)
Vektor
•
•
•
tenzor 1-ého rádu (stupňa)
pre úplné určenie veličiny je potrebné poznať veľkosť, smer a nositeľku
(voľný, viazaný, polohový, jednotkový)
Skalár
•
•
tenzor 0-tého rádu (stupňa)
je určený veľkosťou
prednáška č.2 - 3/59
Základná terminológia
• pri uvažovaní 3D súradnicového systému je tenzorom 2-ého rádu
napr. napätie v bode telesa
• tenzorom 1-ého rádu je napr. sila
• materiálové vlastnosti, ...
prednáška č.2 - 4/59
Základné matematické operácie
Systém lineárnych algebraických rovníc
a11 x1 + a12 x2 + a13 x3 + ... + a1n xn = b1
a21 x1 + a22 x2 + a23 x3 + ... + a2n xn = b2
am1 x1 + am2 x2 + am3 x3 + ... + amn xn = bm
kde x1, x2, ... , xn = sú neznáme,
môžeme v maticovom tvare zapísať
Ax=b
prednáška č.2 - 5/59
Základné matematické operácie
kde
(m  n)
 a11
a
A  [aij ]   21
 

 am1
a12
a22

am 2
 a1n 
 a2 n 

  

 amn 
( n  1)
( m  1)
 x1 
x 
x  [ xi ]   2 

 
 xn 
 b1 
b 
b  [bi ]   2 

 
 xm 
aij - prvok matice na i-tom riadku a j-tom stĺpci
Ak:
m = n - štvorcová matica
m = 1 - riadková matica (vektor)
n = 1 - stĺpcová matica (vektor)
prednáška č.2 - 6/59
Základné matematické operácie
Sčítavanie (odčítavanie) matíc:
! matice A, B musia mať rovnaký rozmer (m  n)
C=A+B
C=A-B
[cij] = [aij] + [bij]
[cij] = [aij] – [bij]
Násobenie matíc:
kA = C
(ln )
[k aij] = [cij]
(lm ) ( mn )
C  A
B
m
cij   aik bkj  aik bkj
i = 1, 2, ... l
asociatívnosť:
avšak:
A (B C) = (A B) C
AB  BA
k 1
j = 1, 2, ... n
prednáška č.2 - 7/59
Základné matematické operácie
Transponovaná matica:
ak:
A = [aij]
potom:
AT = [aji]
Symetrická matica:
štvorcová (n  n) matica A sa nazýva symetrická ak platí:
A = AT [aij] = [aji]
Jednotková matica:
IA = AI = A
( n  n)
1
0
I


0
0  0
1  0

  

0  1
prednáška č.2 - 8/59
Základné matematické operácie
Determinant matice:
Cramerovo pravidlo
n
det A   (1)1 j a1 j det A1 j
j 1
( 22)
a b 
det A  A  det 
 ad  bc

c d 
(33)
a1
det A  det b1

 c1
a2
b2
c2
a3 
b3   a1b2 c3  a2b3c1  a3b1c2  a3b2 c1  a1b3c2  a2b1c3

c3 
Singulárna matica:
det A = 0
prednáška č.2 - 9/59
Základné matematické operácie
Pozitívne definitná matica:
štvorcová (n  n) matica A sa nazýva pozitívne definitná ak pre
ľubovolný nenulový vektor x platí:
xTA x > 0
Matica A je potom nesingulárna.
Inverzná matica:
! existuje iba pre štvorcové a nesingulárne matice, ak det A ≠ 0
A A-1 = A-1A = I
A 1 
1
C T kde
det A
C = [cij]
cij = (-1)i+j Mij
Mij – je determinant matice, ktorá vznikne vynechaním i-teho riadku
a j-teho stĺpca z matice A.
prednáška č.2 - 10/59
Základné matematické operácie
Derivovanie a integrovanie matíc:
d
 daij 
A

dt
 dt 
 Adt   a dt 
ij
Okrem toho platí:
(A B)T = BTAT
(A B C)T = CT(A B)T = CTBTAT
pre štvorcové matice:
A B = BTA
(A B)-1 = B-1A-1
 ( x T A x)
 2A x
x
 (x T b)
b
x
prednáška č.2 - 11/59
Riešenie systému lineárnych
algebraických rovníc
Na vyriešenie systému lineárnych algebraických rovníc môžeme
použiť
Priame metódy
• inverznú maticu
• Gaussovu eliminačnú metódu (Gauss Elimination)
• Gauss-Jordanovu metódu (Gauss-Jordan Elimination)
• Cramerove pravidlo (Cramer´s Rule)
• LU rozklad (LU Decomposition)
• Choleskyho rozklad (Cholesky Decomposition)
Nepriame metódy
• iteračné metódy (približné) – napr. Gauss-Seidelovu metódu,
Jacobiho metódu – vhodné pre riešenie veľkých sústav rovníc
• ...
prednáška č.2 - 12/59
Príklad
Nájdite riešenie (x1 až x3) sústavy rovníc
–x1 + 3x2 – 2x3 = 2
2x1 – 4x2 + 2x3 = 1
4x2 – x3 = 3
Sústavu je možné v maticovom tvare zapísať
 1 3  2  x1  2
 2  4 2   x   1 

 2   
 0
4  1  x3  3
(33) (31)
A
(31)
x  b
prednáška č.2 - 13/59
Inverzná matica:
Ax=b
A-1A x = A-1b
I x = A-1b
x = A-1b
Matica kofaktorov:

4
11
(

1
)

4


3
2 1
C   (1)
4


3
3 1
(1)
4

2
1
2
1
2
2
  4  5  2
CT   2
1  2


 8
4  2
1
CT
det A
C  [(-1)i  j M ij ]
A 1 
2
(1)
0
1
2 2
(1)
0
1
3 2
(1)
2
1 2
2
1
2
1
2
2
2 4 
(1)

0 4 
8
 4 2


1
3

(1) 2  3
4
   5 1

0 4 
 2  2  2


1
3
(1) 3 3

2  4 
1 3
prednáška č.2 - 14/59
Inverzná matica:
det A = |A| = -6
0,8 3 0, 3 
 4  5  2  0, 6
1
1 
A 1  CT 
2
1  2   0, 3  0,1 6 0, 3 

A
6 
 8
4  2   1, 3  0, 6 0, 3 
Výpočet koreňov:
x = A-1b
0,8 3 0, 3  2  0, 6 .2  0, 3 .1  1, 3 .3   3,1 6 
 x1   0, 6
 x    0, 3  0,1 6 0, 3  1  1,8 3 .2  0,1 6 .1  0, 6 .3   0,1 6 
  
 

 2 
 x3    1, 3  0, 6 0, 3  3  0, 3 .2  0, 3 .1  0, 3 .3   2, 3 
prednáška č.2 - 15/59
Inverzná matica:
Vzhľadom na náročný výpočet inverznej matice, z dôvodu nutnosti počítať
determinanty matíc, sa táto metóda v počítačovej mechanike prakticky
nepoužíva.
(viď Cramerova metóda)
prednáška č.2 - 16/59
Cramerovo pravidlo:
( nn ) ( n1)
A
( n1)
n
a x  b
x  b
j 1
ij
j
i
D(i )
xi 
A
kde D(i) je matica, ktorej i-ty stĺpec je nahradený vektorom b.
Výpočet koreňov:
2 3 2
1 4 2
D (1)
3 4 1
 19
x1 


 3,1 6

1
3

2
A
6
2 4 2
0
4
1
prednáška č.2 - 17/59
Cramerovo pravidlo:
1 2  2
2 1 2
D( 2)
0 3 1
1
x2 


 0,1 6
1 3  2  6
A
2 4 2
0
4
1
1 3 2
2 4 1
D( 3)
0
4 3
14
x3 


 2, 3

1
3

2
A
6
2 4 2
0
4
1
prednáška č.2 - 18/59
Cramerovo pravidlo:
Praktické použitie tejto metódy je obmedzené na malé matice n  m (približne
n,m ≤ 5).
Napr. pre výpočet determinantu matice 10x10 by bolo potrebné vykonať cez
30 miliónov operácií. Determinant by mal 10! = 3 628 800 členov.
Preto sa táto metóda v počítačovej mechanike, podobne ako metóda inverznej
matice, nepoužíva.
prednáška č.2 - 19/59
Gaussova eliminačná metóda:
 a11
a
 21
a31
a12
a22
a32
a13   x1  b1 
a23   x2   b2 
   
a33   x3  b3 
 a11
a
 21
a31
a12
a22
a32
a13 b1 
a23 b2   [ A : b]

a33 b3 
tzv.
rozšírená
matica
 1 3  2  x1  2
 2  4 2   x   1 

 2   
 0
4  1  x3  3
prednáška č.2 - 20/59
Gaussova eliminačná metóda:
1. krok
a11

  a a  a
 a 11 21
 aa a11  a31
21
11
31
11
 1
  2 (1)  2
 1
 01 (1)  0
a11
0

 0
a12
a221
a321
a12
 a
 a
 a22
 a 21
a11
12
 a 31
a11
 a32
12


b1
  x1  

a13  a23   x2    aa b1  b2 
a13  a33   x3   aa b1  b3 
a13
 a 21
a11
 a 31
a11
21
11
31
11
2
  x1   2 
  x     2 2  1
2
2
1 3  4
1 ( 2)  2
  2   1

0
0
0
1 3  4
1 ( 2)  1 
  x3   1 2  3
3
a13   x1   b1 
a231   x2   b21 
   
1
a33   x3  b31 
 1 3  2  x1  2
 0 2  2  x   5 

 2   
 0 4  1  x3  3
prednáška č.2 - 21/59
Gaussova eliminačná metóda:
 a11
a
 21
a31
a12
a22
a32
a13   x1  b1 
a23   x2   b2 
   
a33   x3  b3 
a11

  a a  a
 a 11 21
 aa a11  a31
21
11
31
11
a12
 a
 a
 a 21
a11
 a22
12
 a 31
a11
 a32
12
úprava druhého riadku matice


b1
  x1  

a13  a23   x2    aa b1  b2 
a13  a33   x3   aa b1  b3 
a13
 a 21
a11
 a 31
a11
21
11
31
11
prednáška č.2 - 22/59
Gaussova eliminačná metóda:
 a11
a
 21
a31
a12
a22
a32
a13   x1  b1 
a23   x2   b2 
   
a33   x3  b3 
a11

  a a  a
 a 11 21
 aa a11  a31
21
11
31
11
a12
 a
 a
 a 21
a11
 a22
12
 a 31
a11
 a32
12


b1
  x1  

a13  a23   x2    aa b1  b2 
a13  a33   x3   aa b1  b3 
a13
 a 21
a11
 a 31
a11
21
11
31
11
prednáška č.2 - 23/59
Gaussova eliminačná metóda:
 a11
a
 21
a31
a12
a22
a32
a13   x1  b1 
a23   x2   b2 
   
a33   x3  b3 
a11

  a a  a
 a 11 21
 aa a11  a31
21
11
31
11
a12
 a
 a
 a 21
a11
 a22
12
 a 31
a11
 a32
12


b1
  x1  

a13  a23   x2    aa b1  b2 
a13  a33   x3   aa b1  b3 
a13
 a 21
a11
 a 31
a11
21
11
31
11
prednáška č.2 - 24/59
Gaussova eliminačná metóda:
 a11
a
 21
a31
a12
a22
a32
a13   x1  b1 
a23   x2   b2 
   
a33   x3  b3 
a11

  a a  a
 a 11 21
 aa a11  a31
21
11
31
11
a12
 a
 a
 a 21
a11
 a22
12
 a 31
a11
 a32
12


b1
  x1  

a13  a23   x2    aa b1  b2 
a13  a33   x3   aa b1  b3 
a13
 a 21
a11
 a 31
a11
21
11
31
11
prednáška č.2 - 25/59
Gaussova eliminačná metóda:
1. krok
a11

  a a  a
 a 11 21
 aa a11  a31
21
11
31
11
 1
  2 (1)  2
 1
 01 (1)  0
a11
0

 0
a12
a221
a321
a12
 a
 a
 a22
 a 21
a11
12
 a 31
a11
 a32
12


b1
  x1  

a13  a23   x2    aa b1  b2 
a13  a33   x3   aa b1  b3 
a13
 a 21
a11
 a 31
a11
21
11
31
11
2
  x1   2 
  x     2 2  1
2
2
1 3  4
1 ( 2)  2
  2   1

0
0
0
1 3  4
1 ( 2)  1 
  x3   1 2  3
3
a13   x1   b1 
a231   x2   b21 
   
1
a33   x3  b31 
 1 3  2  x1  2
 0 2  2  x   5 

 2   
 0 4  1  x3  3
prednáška č.2 - 26/59
Gaussova eliminačná metóda:
2. krok
a11

0
 0

 1
0

 0
  x1   2 
2
 2   x2    5 
  

4
4
4
2 2 4
2 ( 2)  1
  x3   2 5  3
a11
0

 0
a12
 a 132
a 122
a12
a221
a221  a321
3
a221
0

 a 132
a 122
  x1  

a13
b1
  

1
a231
x

b
2
 2  

1
1
1
1
a
a23  a33   x3   a b2  b3 
1
32
1
22
2
a13   x1   b1 
a231   x2   b21 
   
2
a33   x3  b32 
 1 3  2  x1   2 
 0 2  2  x    5 

 2   
 0 0 3   x3   7
prednáška č.2 - 27/59
Gaussova eliminačná metóda:
Spätnou substitúciou
bnn 1
xn  n 1
ann
- posledný koreň n
1  i 1 n i 1 
xi   bi   air xr 
aii 
r i 1

- ostatné korene (n-1) až 1
Výpočet koreňov:
b32  7
x3  2 
 2, 3
a33
3


x2 
1 1 1
1
b2  a23 x3  5  (2).(2, 3 )   0,1 6
a22
2
x1 
1
b1  a12 x2  a13 x3   1 2  3.0,1 6  (2).(2, 3)   3,1 6
a11
1
prednáška č.2 - 28/59
Gaussova eliminačná metóda:
Táto metóda v počítačovej mechanike používa veľmi často.
prednáška č.2 - 29/59
LU dekompozícia:
A=LU
A x = L U x = L (U x) = b
Ly= b U x=y
 a11
a
 21
a31
A
a12
a22
a32

L
0
a13 
l11
a23   l21 l22


a33 
l31 l32
0
0

l33 
U
u12
u11
0 u
22

 0 0
u13 
u23 

u33 
l11u12
l11u13
 a11 a12 a13  l11 0 0  u11 u12 u13  l11u11

a
  l l
 0 u
  l u l u  l u

a
a
0
u
l
u

l
u
21
22
23
21
22
22
23
21
11
21
12
22
22
21
13
22
23

 

 

 a31 a32 a33  l31 l32 l33   0 0 u33  l31u11 l31u12  l32u22 l31u13  l32u23  l33u33 
l11  1 l22  1 l33  1
prednáška č.2 - 30/59
LU rozklad:
 a11
a
 21
a31
a12
a22
a32
a13   x1  b1 
a23   x2   b2 
   
a33   x3  b3 
 1 3  2  x1  2
 2  4 2   x   1 

 2   
 0
4  1  x3  3
u12
u13
 a11 a12 a13   1 0 0 u11 u12 u13   u11

a
  l
 0 u
  l u l u  l u

a
a
1
0
u
l
u

l
u
22
23 
21 13
22 23
 21 22 23   21

 21 11 21 12 22 22

a31 a32 a33  l31 l32 1  0 0 u33  l31u11 l31u12  l32u22 l31u13  l32u23  l33u33 
prednáška č.2 - 31/59
LU rozklad:
l11  1 l22  1 l33  1
l11u11  a11 l11u12  a12 l11u13  a13
a11
 a11  1
l11
a
 u12  12  a12  3
l11
a
 u13  13  a13  2
l11
a
2
 l21  21 
 2
u11  1
l11u11  a11  u11 
l11u12  a12
l11u13  a13
l21u11  a21
l21u12  l22 u 22  a22

l21u13  l22 u 23  a23
 u 23  a23  l21u13  2  (2).(2)  2
l31u11  a31
 l31 
u 22  a22  l21u12  4  (2).3  2
a31 0

0
u11  1
l31u12  l32 u 22  a32
l31u13  l32 u 23  l33u33  a33
a32  l31u12 4  0.3

2
u 22
2
a l u l u
 1  0.(2)  2.(2)
 u33  33 31 13 32 23 
3
l33
1 prednáška č.2 - 32/59
 l32 
LU rozklad:
Ly =b
l11 0
l l
 21 22
l31 l32
0   y1  b1 
0   y2   b2 
   
l33   y3  b3 
1 0
l
1
 21
l31 l32
0  y1  b1 
0  y2   b2 
   
1  y3  b3 
 1 0 0  y1  2
  2 1 0   y   1 

 2   
 0 2 1  y3  3
Doprednou substitúciou:
 y1   2 
y    5 
 2  
 y3   7
prednáška č.2 - 33/59
LU rozklad:
Ux=y
u11 u12
0 u
22

 0 0
u13   x1   y1 
u23   x2    y2 
   
u33   x3   y3 
 1 3  2  x1   2 
 0 2  2  x    5 

 2   
 0 0 3   x3   7
Spätnou substitúciou:
 x1   3,1 6 
 x    0,1 6 

 2 
 x3   2, 3 
prednáška č.2 - 34/59
LU rozklad:
Skúška správnosti:
Ax=b
 1 3  2  3,1 6  2
 2  4 2   0,1 6   1
  


 0
4  1  2, 3  3
prednáška č.2 - 35/59
LU rozklad:
Táto metóda v počítačovej mechanike používa veľmi často.
prednáška č.2 - 36/59
Choleskyho dekompozícia:
metóda použiteľná len pre symetrické matice
matica musí byť „dostatočne“ pozitívne definitná
A = UT U
 a11
a
 21
a31
A

UT
a12 a13 
0
u11 0
a22 a23   u12 u22 0 
u13 u23 u33 
a32 a33 
u11
0

 0
U
u12 u13 
u22 u23 
0 u33 
2

u11u12
u11u13
 a11 a12 a13   u11

a a
  u u
2
2
a
u

u
u
u

u
u
12
22
12 13
22 23 
 21 22 23   12 11
2
2 
a31 a32 a33  u13u11 u13u12  u23u22 u132  u23
 u33

alternatívne:
A = L LT
prednáška č.2 - 37/59
Choleskyho rozklad:
 6  4 0  x1  1
  4 3 1   x    2

 2   
 0
1 6  x3  3
 a11 a12 a13   x1   b1 
a
  x   b 
a
a
 21 22 23   2   2 
 a31 a32 a33   x3  b3 
Ax=b
prednáška č.2 - 38/59
Choleskyho rozklad:
2
u11
 a11
u11u12  a12
u11u13  a13
 u11  a11  6  2,4495
a
4
 u12  12 
 1,6330
u11 2,4495
a
0
 u13  13 
0
u11 2,4495
2
2
u12
 u 22
 a22
u12 u13  u 22 u 23  a23
2
2
2
u13
 u 23
 u33
 a33
 u 22  a22  u12  3  (1,633) 2  0,5773
a u u
1  (1,633).0
 u 23  23 12 13 
 1,73205
u 22
0,5773
2
2
2
 u33  a33  u13
 u 23
 6  0 2  (0,9428) 2  1,73205
prednáška č.2 - 39/59
Choleskyho rozklad:
u11 u12
U   0 u 22
 0
0
u13 
u 23 
u33 
0 
2,4495  1,6330
U   0
0,5773 1,73205
 0
0
1,73205
 a11
a
 21
a31
A

UT
a12 a13 
0
u11 0
a22 a23   u12 u22 0 
u13 u23 u33 
a32 a33 
u11
0

 0
U
u12 u13 
u22 u23 
0 u33 
A

UT
0
0 
 6  4 0
 2,4495
 4 3 1   1,6330 0,5773

0




 0
 0
1 6
1,73205 1,73205
U
0 
2,4495  1,6330
 0

0
,
5773
1
,
73205


 0
0
1,73205
prednáška č.2 - 40/59
Choleskyho rozklad:
u11
u
 12
u13
0
u22
u23
0   y1   b1 
0   y2   b2 
u33   y3  b3 
0
0   y1  1
 2,4495
 1,6330 0,5773
  y    2
0

 2   
 0
1,73205 1,73205  y3  3
Doprednou substitúciou:
 y1   0,40825 
 y    4,6188 
 2 

 y3   2,88675
prednáška č.2 - 41/59
Choleskyho rozklad:
u11 u12
0 u
22

 0
0
u13   x1   y1 
u23   x2    y2 
u33   x3   y3 
0   x1   0,40825 
2,4495  1,6330
 0
  x    4,61921 
0
,
5773
1
,
73205

 2  

 0
0
1,73205  x3   1,61921
Spätnou substitúciou:
 x1   8,83 
 x    13 

 2 
 x3   1, 6


prednáška č.2 - 42/59
Choleskyho rozklad:
Skúška správnosti:
Ax=b
 6  4 0  8,83  1 
 4 3 1  13   2
  


 0
1 6  1, 6 3
prednáška č.2 - 43/59
Choleskyho rozklad:
Táto metóda v počítačovej mechanike používa veľmi často.
prednáška č.2 - 44/59
Jacobiho a Gauss-Seidelova iteračná metóda:
Majme sústavu rovníc
Ax=b
Pre maticu A (n x n) platí:
aii ≠ 0 pre i = 1...n
Maticu A rozložme
A = AL + AD + AU
AL – dolná trojuholníková matica
AD – diagonálna matica diag[aii]
AU – horná trojuholníková matica
príp. ALT transponovaná dolná trojuholníková matica
Ďalej uvažujeme len prípad:
teda:
AU = ALT
(pre symetrické matice)
A = AL + AD + ALT
prednáška č.2 - 45/59
Jacobiho a Gauss-Seidelova iteračná metóda:
 a11
a
 21
 a31

 
 an1
a12
A
a13
a22
a23
a32

a33

an 2
an 3

 a1n 
0
a
 a1n 
 21
 a1n    a31


  
 
 an1
 ann 
AL
0
0
0
0
a32

0

an 2
an 3

 0
a11
0
 0

 0   0


 
 
 0
 0

0
AD
0

a22
0

0

0
a33 
 
0

0 
0 a12
0 0
0 

0   0 0


 
 
0 0
ann 
A TL
a13  a1n 
a23  a1n 
0  a1n 

   
0  0 
prednáška č.2 - 46/59
Jacobiho a Gauss-Seidelova iteračná metóda:
Iteračný algoritmus pre výpočet koreňov
Jacobiho iteračná metóda
x ( k 1)

1 

bi 
aii 





aij x (jk ) 
j 1

j i

n

x ( k 1)  A D1 b  ( A L  A TL )x ( k )

x ( k 1)  x ( k )
Gauss-Seidelova iteračná metóda
i 1
n

1 
( k 1)
( k 1)
(k )
x
 bi   aij x j   aij x j 
aii 
j 1
j i 1


x ( k 1)  A D1 b  A L x ( k 1)  A TL x ( k )
x
(k )


prednáška č.2 - 47/59
Jacobiho a Gauss-Seidelova iteračná metóda:
Iteračný algoritmus pre výpočet koreňov
Jacobiho iteračná metóda
M J x ( k 1)  N J x ( k )  b
kde
MJ  AD
N J  ( A L  A TL )
Gauss-Seidelova iteračná metóda
M G x ( k 1)  N G x ( k )  b
kde
MG  AL  AD
N G   A TL
prednáška č.2 - 48/59
Jacobiho a Gauss-Seidelova iteračná metóda:
Príklad:
Jacobiho metódou nájdite korene sústavy rovníc s presnosťou
 ≤ 0,0003
Ax=b
 6  4 0  x1  1 
 4 3 1   x    2

 2   
 0
1 6  x3  3
Rozložíme
maticu
A
 A AL
6
 4

 0
 a11
a
 21
a31
4
3
1
a12
a22
a32
0
0
1   4
 0
6
a13 
0
a23   a21
a31
a33 
0
0
1
0
0
a32

D

0
6 0 0 
0
0 3 0 
0
0 




0 0 6
0
0
0
0
a11 0
0
0   0 a22 0   0
 0
0
0
0 a33 
A TL
4
0
0
a12
0
0
0
1
0
a13 
a23 
0 
prednáška č.2 - 49/59
Jacobiho a Gauss-Seidelova iteračná metóda:
x(k 1)  MJ1 NJ x(k )  MJ1 b
6 0 0 
M J  0 3 0
0 0 6

  0 0 0  0

N J    4 0 0  0
  0 1 0  0
 

 16

M J1  0
0

 4 0  0
 

0 1   4
0 0  0
0 0

1
0
3

0 16 
0
 1
 1 0 
4
0
 16 0 0 0 4 0  0 23



M J1 N J  0 13 0 4 0  1   43 0
0 0 16  0  1 0  0  16



0

 13 
0 
 16 0 0 1   16 


 
M J1 b  0 13 0 2   23 
0 0 16  3  12 
 


prednáška č.2 - 50/59
Jacobiho a Gauss-Seidelova iteračná metóda:
x ( k 1)
0 23

  43 0
0  1
6

 16 
0

 
 13  x ( k )   23 
 12 
0 
 
Štartovacie riešenie (zvolené riešenie)
x (0)
0 
 0
0
1. iterácia
x (1)
0 23

  43 0
0  1
6

 16  0 23
0

  
 13  x ( 0)   23    43 0
 12  0  16
0 
  
0  0  16   16  0,16

2 2 
1  
 3  0   3    3    0, 6 
0  0  12   12   0,5 
prednáška č.2 - 51/59
Jacobiho a Gauss-Seidelova iteračná metóda:
2. iterácia
x ( 2)
0

  43
0

x ( 2)  x (1)
x
(1)
 23
0
1
6

 16   0
0
2  4
1  (1)
x

3
 3    3
 12   0
0
  
0
1
6
11
  0,61 
0  16   16   18

 2   13  
1 2
  3    18   0,7 2
3 3
0  12   12   187   0,38 
( x1( 2)  x1(1) ) 2  ( x2( 2)  x2(1) ) 2  ( x3( 2)  x3(1) ) 2
( x1(1) ) 2
 ( x1(1) ) 2
11
11
( 185  16 ) 2  ( 18
 23 ) 2  ( 18
 12 ) 2
( 16 ) 2
 23
11 2
 ( 18
)
11 2
 ( 18
)
 ( x1(1) ) 2

 0,543021
prednáška č.2 - 52/59
Jacobiho a Gauss-Seidelova iteračná metóda:
3. iterácia
x ( 3)
0

  43
0

x ( 3)  x ( 2 )
x
( 2)
 23
0
1
6

 16   0
0
2  4
1  ( 2)
x

3
 3    3
 12   0
0
  
 23
0
1
6
35
11
  16   54
 0,648148
0  18
 2   73  

1   13 



1
,
35185
3   18 
 3   54  

41 
0  187   12   108
  0,37963 
( x1(3)  x1( 2) ) 2  ( x2(3)  x2( 2) ) 2  ( x3(3)  x3( 2) ) 2
( x1( 2) ) 2
 ( x1( 2) ) 2
 ( x1( 2) ) 2
 0,616673  
prednáška č.2 - 53/59
Jacobiho a Gauss-Seidelova iteračná metóda:
4. iterácia
x ( 4)
0

  43
0

x ( 4 )  x ( 3)
x
( 3)
 23
0
1
6

 16   0
0
2  4
1  ( 3)
x

3
 3    3
 12   0
0
  
 23
0
1
6
35
173
  16   162
  1,0679 
0  54
 2   455  

1   73 



1
,
40432
3   54 
 3   324  

89
41
1
0  108   2   324  0,274691
( x1( 4)  x1(3) ) 2  ( x2( 4)  x2(3) ) 2  ( x3( 4)  x3(3) ) 2
( x1(3) ) 2
 ( x1(3) ) 2
 ( x1(3) ) 2
 0,281821  
prednáška č.2 - 54/59
Jacobiho a Gauss-Seidelova iteračná metóda:
164. iterácia
x (164 )
0
 4
  3

0
x (164 )  x (163)
x
(163 )
 23
0
1
6

 16   8,75311 
0
 (163)  2  

1
x
  3    12,8802 
3

1 

0

1
,
64661


2
( x1(164 )  x1(163) ) 2  ( x2(164 )  x2(163) ) 2  ( x3(164 )  x3(163) ) 2
( x1(163) ) 2  ( x1(163) ) 2  ( x1(163) ) 2
 0,000287561 
Presné riešenie:
 x1   8,83 
 x    13 

 2 
 x3   1, 6


prednáška č.2 - 55/59
Jacobiho a Gauss-Seidelova iteračná metóda:
Skúška správnosti:
Ax=b
 6  4 0  8,75311  0,997832


 


4
3
1
12
,
8802

1
,
98157


 

0
1 6  1,64661  3,00054 

Presné riešenie:
 6  4 0  8,83  1 
 4 3 1  13   2
  


 0
1 6  1, 6 3
prednáška č.2 - 56/59
Zle podmienené matice:
Majme sústavu rovníc
Ax=b
Ak je matica A zle podmienená (ill-conditioned matrix), riešenie úlohy
môže byť obtiažne. Napr. malé zaokrúhlenia členov b vektora, môže
výrazne ovplyvniť riešenie /korene/ sústavy rovníc (x vektor).
prednáška č.2 - 57/59
Zle podmienené matice:
Príklad:
Riešením sústavy rovníc
9x + 8y = 0,8
8x + 7y = 0,7
sú korene:
x=0
y = 0,1
Ak zavedieme malú chybu do pravej strany rovníc (vektor b) (napr.
zaokrúhlovacou chybou)
9x + 8y = 0,81
8x + 7y = 0,69
sú korene:
x = -0,15
y = 0,27
prednáška č.2 - 58/59
Riešenie systému nelineárnych
algebraických rovníc
prednáška č.2 - 59/59
Riešiče implementované
®
v programe ΛNSYS
Direct Solvers - priame riešiče
• Sparse Direct Solver – využíva LU rozklad
• Frontal (Wavefront) Solver – založený na metóde Gaussovej eliminácie
Iterative Solver - iteratívne riešiče
•
•
•
•
Jacobi Conjugate Gradient (JCG) Solver
Preconditioned Conjugate Gradient (PCG) Solver
Incomplete Cholesky Conjugate Gradient (ICCG) Solver
Automatic Iterative (Fast) Solver Option
Parallel/Distributed Solvers - distribuované riešiče
• Algebraic Multigrid (AMG) Solver
• Distributed Jacobi Conjugate Gradient (DJCG) Solver
• Distributed Preconditioned Conjugate Gradient (DPCG) Solver
prednáška č.2 - 60/59