Metody rozwiązywania układów równań

liniowych

n

n

nn

n

n

n

n

n

n

b

x

a

x

a

x

a

b

x

a

x

a

x

a

b

x

a

x

a

x

a





























2

2

1

1

2

2

2

22

1

21

1

1

2

12

1

11

b

Ax





















































n

nn

n

n

n

n

b

b

b

a

a

a

a

a

a

a

a

a

















2

1

2

1

2

22

21

1

12

11

b

A

0

det 

A

Parę słów o macierzach

Macierz m



n: tablica m na n (m wierszy n kolumn) liczb (np.

tabliczka mnożenia).

Macierz kwadradowa: m=n

Macierz symetryczna (zawsze kwadratowa): a

ij

=a

ji

Macierz transponowana A

T

: (A

T

)

ij

=a

ji

Macierz nieosobliwa: macierz o niezerowym wyznaczniku.

Macierz dodatnio określona:

x

T

Ax>0 dla każdego niezerowego wektora x.

Norma euklidesowa macierzy:

Norma spektralna macierzy

Wskaźnik uwarunkowania macierzy









n

i

n

j

ij

a

1

1

2

A

i

i



max



A

1

cond





A

A

A

Metody skończone:

• Metoda Gaussa
• Metoda Gaussa-Jordana
• Metody Choleskiego
• Metoda Householdera
• Metoda sprzężonych gradientów

Metody iteracyjne dla dużych układów równań:

• Metoda Jacobiego
• Metoda Gaussa-Seidla

n

n

nn

n

n

n

n

c

x

r

c

x

r

x

r

c

x

r

x

r

x

r























2

2

2

22

1

1

2

12

1

11

Metoda eliminacji Gaussa z wyborem elementu głównego w
kolumnie

Układ równań sprowadzamy do postaci trójkątnej

Układ z macierzą trójkątną można następnie łatwo rozwiązać
zaczynając od obliczenia wartości x

n

z n-tego równania,

następnie wstawić x

n

do równania n-1 i wyliczyć z niego x

n-1

,

następnie wstawić x

n

oraz x

n-1

do równania n-2 i wyliczyć x

n-2

aż

do dotarcia do równania pierwszego i wyznaczenia x

1

.

1. Wybieramy równanie i takie, że |a

i1

| jest największym elementem w

pierwszej kolumnie po czym przestawiamy i-te równanie na początek i
eliminujemy x

1

z równań od 2 do n.

)

1

(

)

1

(

2

2

)

1

(

2

)

1

(

2

)

1

(

2

2

)

1

(

22

)

0

(

1

)

0

(

1

2

)

0

(

12

1

)

0

(

11

n

n

n

n

n

n

n

n

b

x

a

x

a

b

x

a

x

a

b

x

a

x

a

x

a































1

1

,

1

)

0

(

11

)

0

(

1

)

0

(

1

)

0

(

)

1

(

)

0

(

11

)

0

(

1

)

0

(

1

)

0

(

)

1

(















i

a

a

b

b

b

k

i

a

a

a

a

a

i

i

i

i

k

ik

ik

2. Procedurę powtarzamy z macierzą A

(1)

o rozmiarach (n-1)x(n-1) i

wektorem b

(1)

o rozmiarze n-1, eliminując z nich drugą zmienną i

otrzymując macierz A

(2)

o rozmiarach (n-2)x(n-2) i wektor b

(2)

o

rozmiarze n-2. W ten sam sposób postępujemy z kolejnymi macierzami
A

(2)

, A

(3)

,..., A

(n-1)

oraz wektorami b

(2)

, b

(3)

,..., b

(n-1)

.

j

i

a

a

b

b

b

j

k

j

i

a

a

a

a

a

j

jj

j

ij

j

j

j

i

j

i

j

jj

j

ij

j

jk

j

ik

j

ik































)

1

(

)

1

(

)

1

(

)

1

(

)

(

)

1

(

)

1

(

)

1

(

)

1

(

)

(

,

Dla j-tego kroku

Po zakończeniu operacji otrzymujemy układ równań z macierzą
trójkątną

)

1

(

)

1

(

22

)

0

(

11

)

1

(

)

1

(

)

2

(

3

)

0

(

3

3

)

2

(

33

)

1

(

2

)

0

(

2

3

)

1

(

23

2

)

1

(

22

)

0

(

1

)

0

(

1

3

)

0

(

13

2

)

0

(

12

1

)

0

(

11

)

1

(

det





































n

nn

p

n

n

n

n

nn

n

n

n

n

n

n

a

a

a

b

x

a

b

x

a

x

a

b

x

a

x

a

x

a

b

x

a

x

a

x

a

x

a















A

p jest liczbą przestawień wierszy macierzy A podczas
sprowadzania układu równań do postaci trójkątnej.

1

,...,

2

,

1

1

)

1

(

)

1

(

)

1

(

)

1

(

)

1

(

)

1

(































n

n

j

x

a

a

a

b

x

a

b

x

k

n

j

k

j

jj

j

jk

j

jj

j

j

j

n

nn

n

n

n

3. Z otrzymanego układu równań z macierzą trójkątną

wyznaczamy po kolei x

n

, x

n-1

,..., x

1

.

Wysiłek obliczeniowy (liczba mnożeń i dzieleń) w metodzie
eliminacji Gaussa:

Faktoryzacja macierzy A: n(n

2

-1)/3 operacji

Przekształcenie wektora b: n(n-1)/2 operacji

Obliczenie x: n(n+1)/2 operacji.

Razem: n

3

/3+n

2

-n/3≈n

3

/3 operacji.

Kod źródłowy metody eliminacji Gaussa.

Zastosowania metody eliminacji Gaussa i innych metod
rozwiązywania układów równań liniowych.

1. Odwracanie macierzy

 

 

 

 

 













































































































0

0

1

0

0

1

,

,

1

,

1









in

i

i

i

i

i

i

i

1

1

1

1

1

1

A

A

A

A

A

A

I

AA

Należy więc jednocześnie rozwiązać n układów równań z wektorami
jednostkowymi poszczególnych osi w przestrzeni n-wymiarowej jako
wyrazami wolnymi. Wszystkie wektory wyrazów wolnych
przekształca się jednocześnie w czasie faktoryzacji macierzy A.
Koszt obliczeniowy operacji wynosi (4/3)n

3

-n/3 operacji.

2. Mnożenie macierzy przez macierz odwrotną

B

AC

B

A

C

1







Jeżeli macierz A ma rozmiar nxn to macierz B musi mieć n
wierszy (liczba kolumn może być dowolna, np. m).

Zagadnienie sprowadza się do jednoczesnego rozwiązania m
układów równań liniowych, których wektory wyrazów wolnych
są równe kolejnym kolumnom macierzy B. Podobie jak przy
odwracaniu macierzy podczas faktoryzacji A jednocześnie
przekształcamy wszystkie kolumny macierzy B.

Wysiłek obliczeniowy wynosi n

3

/3-n/3+mn

2

operacji.

Metoda Gaussa-Jordana

W odróżnieniu od metody eliminacji Gaussa eliminujemy
poszczególne zmienne nie tylko z równań następujących po danym
równaniu ale również z równań je poprzedzających. W wyniku tej
operacji końcowy układ równań ma postać diagonalną a w kroku j-
1 ma on następującą postać:

)

(

)

(

1

)

(

1

,

)

(

2

)

(

1

)

(

1

,

2

)

(

)

(

1

)

(

1

1

)

(

1

,

1

1

)

0

(

11

j

n

n

j

nn

j

j

j

n

j

n

j

jn

j

j

j

j

j

jj

j

n

j

n

j

j

j

b

x

a

x

a

b

x

a

x

a

x

a

b

x

a

x

a

x

a

















































)

1

(

)

1

(

)

1

(

)

1

(

)

1

(

)

1

(

)

(

)

1

(

)

1

(

)

1

(

)

1

(

)

(

,...,

1

,

1

,...,

2

,

1

,...,

1

,

1

,...,

2

,

1

,...,

2

,

1

0

























































i

ii

n

i

i

j

jj

j

ij

j

j

j

i

j

i

j

jj

j

ij

j

jk

j

ik

j

ik

a

b

x

n

j

j

i

a

a

b

b

b

j

i

j

k

a

a

a

a

n

j

j

i

j

k

a

Dla dużych n metoda Gaussa-Jordana wymaga ok. n

3

/2

operacji.

Metody typu Choleskiego dla macierzy symetrycznych silnie
nieosobliwych

T

LDL

A 

L

T

L

D

L

T

LL

A 

klasyczna metoda Choleskiego

tylko dla macierzy dodatnio
określonych.

Postępowanie przy rozwiązywaniu układów równań metodą

faktoryzacji Choleskiego.

1. Wyznaczenie faktorów L i D. Układ przyjmuje postać

LDL

T

x=b

2. Obliczenie pomocniczego wektora w.

w=L

-1

b przez rozwiązanie układu równań Lw=b.

Ponieważ L jest macierzą trójkątną dolną układ ten rozwiązuje się
wyliczając kolejno w

1

, w

2

,…, w

n

podobnie jak w koncowym etapie

eliminacji Gaussa.

3. Obliczenie z=D

-1

w (D jest macierzą diagonalną więc po prostu

dzielimy w

i

przez d

ii

. Ten etap nie występuje w klasycznej metodzie

Choleskiego.

4. Obliczenie x poprzez rozwiązanie układu równań z macierzą

trójkątną górną

L

T

x=z

Ten etap jest identyczny z ostatnim etapem metody eliminacji
Gaussa.

Metoda wymaga ok. n

3

/6 operacji (2 razy mniej niż metoda eliminacji

Gaussa). Uwaga: klasyczna metoda Choleskiego wymaga ponadto
n pierwiastkowań.

Klasyczna faktoryzacja Choleskiego (A=LL

T

)





































1

1

1

1

2

11

1

1

11

11

1

i

k

jk

ik

ij

ii

ji

i

k

jk

ii

ii

j

j

l

l

a

l

l

l

a

l

l

a

l

a

l

Faktoryzacja “bezpierwiastkowa”

kod źródłowy



























































1

0

1

0

0

1

1

,

1

21

2

1

n

n

n

n

l

l

l

d

d

d

















L

D

i

n

k

ik

jk

k

ji

ji

j

j

n

k

nk

k

nn

n

i

k

ki

k

ii

i

d

l

l

d

a

l

d

a

l

l

d

a

d

l

d

a

d

a

d















































1

1

1

1

1

1

1

2

1

1

2

11

1

/

Metoda Householdera

 



























































































)

2

(

)

2

(

1

)

2

(

12

)

2

(

11

)

1

(

)

2

(

)

1

(

1

)

1

(

21

)

1

(

1

)

1

(

11

)

1

(

11

1

1

1

2

1

)

1

(

)

1

(

)

2

(

)

1

(

~

sgn

2

,

A

A

H

A

a

v

v

v

v

I

H

c

b

Q

Ax

Q

cRz

b

Q

Ax

Q

Rz

Q

Q

H

H

H

Q

QR

A

T

T

T

T

T

1

n

n

T

n

a

a

a

a

a

a

a







 





































































































0

~

0

sgn

2

)

1

(

)

1

(

)

1

(

)

2

(

2

)

2

(

22

)

2

(

1

)

2

(

12

)

2

(

11

)

(

)

1

(

)

1

(

)

(

)

(

,

1

)

(

)

(

)

(

2

1

)

(

i

n

in

i

ii

n

n

i

i

i

i

ni

i

i

i

i

i

i

ii

i

ii

i

T

i

i

i

i

n

i

a

a

a

a

a

a

a

a

a

a

a

A

A

H

A

a

v

v

v

v

I

H













Iteracyjna poprawa rozwiązania układu równań liniowych

)

1

(

)

(

)

1

(

)

(

)

(

)

1

(

)

1

(

)

0

(

)

1

(

)

0

(

)

0

(

)

1

(



























p

p

p

p

p

p

z

x

x

r

Ax

b

Az

z

x

x

r

Ax

b

Az

Proces kończymy gdy







)

(

)

1

(

max

max

p

i

i

p

i

i

x

z

Metoda sprzężonych gradientów

Rozwiązanie układu równań przedstawiamy jako rozwiązanie
zagadnienia minimalizacji funkcji E(x), której gradient jest różnicą
pomiędzy lewą a prawą stroną układu równań.









b

Ax

x

b

Ax

A

b

Ax

x

1

T













)

(

grad

2

1

)

(

E

E

Metoda polega na wykonaniu n kroków w kierunkach d

0

, d

1

,…, d

n-1

stopniowo zbliżających się do kierunku jednej z osi
hiperparaboloidy E(x); długość kroku w każdym kierunku jest tak
dobrana aby zlokalizować minimum E na tym właśnie kierunku.
Poszczególne kierunki są A-ortogonalne, tj.

d

k+1

Ad

k

=0

1. Wybieramy punkt startowy x

0

. Obliczamy d

0

=-g

0

=b-Ax

0

. Jeżeli

g

0

=0 procedura jest zakończona.

2. W kroku k=0,1,…,n-1 obliczamy

k

T
k

k

T
k

Ad

d

g

d





k



Wstawiamy x

k+1

=x

k

+

k

d

k

Wstawiamy g

k+1

=g

k

+

k

Ad

k

. Jeżeli ||g

k+|

||< kończymy proces,

w przeciwnym wypadku wstawiamy

k

T
k

k

T

1

k

Ad

d

Ad

g





k



a następnie wstawiamy d

k+1

=-g

k+1

+

k

d

k

Procedura jest teoretycznie zawsze zbieżna w n krokach (0,1,

…,n-1). Dla dobrze uwarunkowanego zadania zbieżność jest
szybsza a dla źle uwarunkowanego zadania często występują
oscylacje.

Kod źródłowy metody sprzężonych gradientów.

Metody iteracji Jacobiego oraz iteracji Gaussa-Seidla

Metody te stosuje się dla dużych układów równań, pojawiających

się przykładowo przy rozwiązywaniu równań różniczkowych
cząstkowych lub równań całkowych.

Iteracja Jacobiego

1. Wyznaczamy przybliżenie początkowe x

(0)

2. Kolejne przybliżenia x

(1)

, x

(2)

,…, obliczamy wyznaczając x

i

z i-tego

równania dla i=1,2,…,n a następnie podstawiając po prawej
stronie poprzednie przybliżenie x.

ii

i

n

i

k

k

p

k

ii

ik

p

i

a

b

x

a

a

x















1

)

(

)

1

(

3. Proces kończymy jeżeli ||x

(p+1)

-x

(p)

||<.

Metoda iteracji Gaussa-Seidla.

Metoda ta różni się od metody iteracji Jacobiego tym, że kolejne
przybliżenie i-tej współrzędnej wektora x obliczamy wykorzystując
przybliżenia x

1

,…,x

i-1

wyliczone w aktualnej iteracji oraz przybliżenia

x

i+1

,…,x

n

z poprzedniej iteracji.

ii

i

i

ii

ik

r

ik

ii

ik

l

ik

i

n

i

k

p

k

r

ik

i

k

p

k

l

ik

ii

i

n

i

k

p

k

ii

ik

i

k

p

k

ii

ik

p

i

a

b

c

i

k

a

a

i

k

b

i

k

i

k

a

a

b

c

x

b

x

b

a

b

x

a

a

x

a

a

x















































































0

0

)

(

)

(

1

)

(

)

(

1

1

)

1

(

)

(

1

)

(

1

1

)

1

(

)

1

(

Metoda iteracji Gaussa-Seidla jest zawsze zbieżna dla macierzy
dodatnio określonych.

)

1

(

)

(

)

(

)

1

(

2

)

(

)

(

)

1

(

2

1

)

(

1

1

)

(

1

1

)

1

(

)

(

)

1

(

max

1

1

1

1

2

1

1

2

















































































p

k

p

k

p

x

p

k

k

p

p

p

opt

ii

i

p

i

n

k

p

i

ii

ik

i

k

p

i

ii

ik

p

i

p

i

x

x

x

x

q

q

q

a

b

x

x

a

a

x

a

a

x

x













Kod źródłowy metody Gaussa-Seidla.

Optymalna wartość czynnika
relaksacji; 

1

jest największą

wartością własną macierzy
B=B

(l)

+B

(r)

Wartość
czynnika
relaksacji łatwa
do obliczenia.