x+y+z = 1, x-y-z = 2, 2x+y-z = 2.
x+y+z = 1, -2y-2z = 1, -y-3z = 0.
x+y+z = 1,
-2y-2z = 1,
-2z = -0.5.
Takaisinsijoitusvaihe (back substitution)
Viimeinen yhtälö => z=0.25.
Sijoitetaan keskimmäiseen yhtälöön: -2y-0.5=1 => y=-0.75.
Sijoitetaan ensimmäiseen yhtälöön: x-0.5=1 => x=1.5. Muunnetaan kerroinmatriisi yläkolmiomatriisiksi, jossa lävistäjän alapuolella on pelkkiä nollia.
! Gaussin eliminointi
! eliminointivaihe
do i=1,n-1
do k=i+1,n
c=A(k,i)/A(i,i)
do j=i,n
A(k,j)=A(k,j)-c*A(i,j)
end do
b(k) = b(k)-c*b(i)
end do
end do
Takaisinsijoitusvaihe etenee alhaalta ylös:
! Gaussin eliminointi
! takaisinsijoitusvaihe
do i=n,1,-1
x(i)=b(i)/A(i,i)
do k=i-1,1,-1
b(k)=b(k)-A(k,i)*x(i)
end do
end do
Kerroinmatriisi A korvautuu yläkolmiomatriisilla;
alkuperäiset kertoimet tuhoutuvat.
Sijoistusvaiheessa käsiteltävän rivin alapuolella olevia vektorin b alkioita ei enää tarvita, joten ratkaisu voidaan tallettaa vektoriin b.
! takaisinsijoitusvaihe
! ratkaisu kirjoitetaan vakiovektorin paalle
do i=n,1,-1
b(i)=b(i)/A(i,i)
do k=i-1,1,-1
b(k)=b(k)-A(k,i)*b(i)
end do
end do
Jos yhtälöitä on n kappaletta,
tilaa tarvitaan n2+n muuttujalle.
Eliminointi ei onnistu, jos jakaja on nolla. Järjestellään yhtälöitä niin, että jakajaksi saadaan nollasta poikkeava luku. Yhtälöiden järjestyksen vaihtaminen ei vaikuta ratkaisuun.
Rivien vaihtamista kutsutaan osittaiseksi pivotoinniksi.
Samalla kannattaa etsiä sellainen rivi, josta jakajaksi saadaan mahdollisimman suuri luku.
Osittaista pivotointia käyttämällä ohjelma on:
! Gaussin eliminointi
! eliminointivaihe osittaisella pivotoinnilla
do i=1,n
! etsitaan suurin alkio sarakkeelta i
m=i
s=A(i,i)
do k=i+1,n
if (abs(A(k,i)) > s) then
s=abs(A(k,i)); m=k
end if
end do
! suurin alkio on rivilla m; vaihdetaan rivit i ja m
v=A(i,:); A(i,:)=A(m,:); A(m,:)=v
x=b(i); b(i)=b(m); b(m)=x
! eliminoidaan muuttuja i
if (A(i,i)==0) exit
do k=i+1,n
c=A(k,i)/A(i,i)
do j=i,n
A(k,j)=A(k,j)-c*A(i,j)
end do
b(k) = b(k)-c*b(i)
end do
end do
Täydellisessä pivotoinnissa jakaja maksimoidaan vaihtamalla lisäksi
pystyrivejä. Tuntemattomien järjestys muuttuu, joten samalla on
pidettävä yllä listaa, jonka avulla ratkaisu osataan lopuksi järjestää
alkuperäistä yhtälöryhmää vastaavaan järjestykseen.
Käytännössä matriisin rivejä ei kannata vaihtaa, vaan käytetään indeksivektoria, joka kertoo kulloisenkin rivien järjestyksen.
! Gaussin eliminointi
! alustetaan indeksivektori
do i=1,n
ind(i)=i
end do
! eliminointivaihe osittaisella pivotoinnilla
do i=1,n
! etsitaan suurin alkio sarakkeelta i
m=i
s=A(ind(m),i)
do k=i+1,n
if (abs(A(ind(k),i)) > s) then
s=abs(A(ind(k),i)); m=k
end if
end do
! suurin alkio on rivilla m;
! vaihdetaan rivien i ja m indeksit
l=ind(i); ind(i)=ind(m); ind(m)=l
! eliminoidaan muuttuja i
ii=ind(i)
if (abs(a(ii,i)) < limit) then exit
do k=i+1,n
kk=ind(k)
c=A(kk,i)/A(ii,i)
do j=i,n
A(kk,j)=A(kk,j)-c*A(ii,j)
end do
b(kk) = b(kk)-c*b(ii)
end do
end do
! takaisinsijoitusvaihe
! ratkaisu kirjoitetaan vakiovektorin paalle
do i=n,1,-1
ii=ind(i)
b(ii)=b(ii)/A(ii,i)
do k=i-1,1,-1
kk=ind(kk)
b(kk)=b(kk)-A(kk,i)*b(ii)
end do
end do
! tulostetaan ratkaisu
do i=1,n
write (6,*) b(ind(i))
end do
(3+2(n-i+1))(n-i)=2(n-i)2+5(n-i).
Nämä suoritetaan muuttujan i arvoilla 1, ... , n-1, joten kaikkiaan operaatioita on
N1 = sum_i=1 .. n-1 [ 2(n-i)2+5(n-i) ]
= sum_i=1 .. n-1 [2n2+5n +2i2-(4n+5)i ]
= (n-1)(2n2+5n)+
2 [ sum_i=1 .. n-1 i2 ] -
(4n+5) [ sum_i=1 .. n-1 i ].
Summaa sum ip voidaan arvioida integraalilla:
sum_i=1 .. n ip \approx int_0^n ip di = np+1/(p+1).
Siis sum_i=1 .. n i \propto n2 ja sum_i=1 .. n i2 \propto n3.
N1 \propto n3.
Takaisinsijoitusvaiheessa sisempi silmukka suoritetaan i-1 kertaa ja kullakin kierroksella tehdään kaksi laskutoimitusta. Tämä silmukka ja yksi jakolasku tehdään muuttujan i arvoilla i, ..., n.
N2 = sum_i=1 .. n (1+2(i-1)) = sum_i=1 .. n (2i-1) \propto n2.
Kaikkiaan liukulukuoperaatioiden määrä
N = N1+N2 \propto n3
Jos menetelmän 1 laskutoimitusten määrä on N1=10n3 ja menetelmän 2 N2=1000n2, kumpikin on yhtä tehokas, kun n=100. Tätä suuremmilla tehtävillä menetelmä 2 on pienemmän eksponenttinsa vuoksi tehokkaampi kuin menetelmä 1. Menetelmän 2 etu kasvaa sitä suuremmaksi, mitä suurempia tehtäviä joudutaan ratkomaan. Jos sen sijaan n on lähes aina pienempi kuin 100, kannattaakin käyttää menetelmää 1, sillä se on ilmeisesti hyvin yksinkertainen ja siksi tehokas pienissä ongelmissa.
Jos laskutoimitusten määrälle voidaan esittää yläraja jonkin tehtävän kokoa kuvaavan parametrin polynomina, tehtävä on ratkaistavissa polynomisessa ajassa.
Jos ajan tarve kasvaa nopeammin kuin mikään polynomi, tehtävä on NP-täydellinen. Jos esimerkiksi tehokkaimmankin ratkaisualgoritmin ajan tarve kasvaa eksponentiaalisesti, se on NP-täydellinen. Tunnetaan useita ongelmia, joille polynomisessa ajassa selviytyvää algoritmia ei ole löydetty. Esimerkiksi eräät kombinatoriikan tehtävät näyttävät olevan NP-täydellisiä.
Matriisi voidaan esittää usealla tavalla ala- ja yläkolmiomatriisien tulona (LU-hajotelma). Hajotelma voidaan tallettaa alkuperäisen kerroinmatriisin paikalle.
a11 a12 ... a1n
a21 a22 ... a2n
. =
.
.
an1 an2 ... ann
l11 0 ... 0 1 u12 ... u1n
l21 l22 ... 0 0 1 ... u2n
. .
. .
. .
ln1 ln2 ... lnn 0 0 ... 1
Esimerkki. Yritetään etsiä hajotelma:
l11 0 0 1 u12 u13 1 2 2 l21 l22 0 0 1 u23 = 1 1 0 l31 l32 l33 0 0 1 2 -1 1Lasketaan kolmiomatriisien tulo:
l11 l11u12 l11u13
l21 l21u12+l22 l21u13+l22u23
l31 l31u12+l32 l31u13+l32u23+l33
Ensinnäkin matriisin L ensimmäinen
pystyrivi on täsmälleen sama kuin alkuperäisessä
matriisissa:
l11 = a11=1,
l21 = a21=1,
l31 = a31=2.
Tulomatriisin ensimmäisen vaakarivin kaksi muuta alkiota antavat yhtälöt
l11u12 = a12,
l11u13 = a13.
Näissä esiintyvä l11 on jo laskettu edellä, joten voimme ratkaista matriisin U ensimmäisen vaakarivin alkiot
u12 = a12/l11 = 2/1=2,
u13 = a13/l11 = 2/1=2.
Toisesta pystyrivistä ovat vielä jäljellä yhtälöt
l21u12+l22 = a22,
l31u12+l32 = a32.
Näissä taaskin kaikki muu on tunnettua paitsi matriisin L toisen pystyrivin alkiot.
l22 = a22-l21u12=1-1\times 2 = -1,
l32 = a32-l31u12=-1-2\times 2 = -5.
Seuraavaksi lasketaan toinen vaakarivi:
l21u13+l22u23=a23.
Tästä saadaan
u23=(a23-l21u13)/l22=(0-1\times 2)/(-1)=2.
Viimeisenä on vuorossa kolmas pystyrivi:
l31u13+l32u23+l33=a33,
josta
l33=a33-l31u13-l32u23 = 1-2\times2-(-5)\times 2=7.
Etsitty hajotelma on siis
1 0 0 1 2 2 1 -1 0 x 0 1 2 2 -5 7 0 0 1Laskentajärjestys on oleellinen. Hajotelman laskemisessa vuorottelevat matriisin L pystyrivit ja matriisin U vaakarivit.
Yhtälöitä tutkimalla havaitaan, että alkioiden lij ja uij laskemisessa tarvitaan alkuperäisestä matriisista vain alkiota aij. Siksi tämä alkio voidaan ilman haittaa korvata lasketulla arvolla lij tai uij.
Kussakin vaiheessa alkioiden lij ja uij laskemiseen käytetään jo laskettuja L ja U-matriisien alkioita, jotka siis on jo talletettu alkuperäisen matriisin alkioiden tilalle.
Hajotelma voidaan tallettaa muodossa
1 2 2 1 -1 2 2 -5 7
! Matriisin LU-hajotelma
! ensimmainen pystyrivi
do i=1,n
L(i,1)=A(i,1)
end do
! ensimmainen vaakarivi
U(1,1)=1
do j=2,n
U(1,j)=A(1,j)/L(1,1)
end do
! muut rivit
do m=2,n
! ensin L:n pystyrivi
do i=m,n
s=0.0
do k=1,m-1
s=s+L(i,k)*U(k,m)
end do
L(i,m)=A(i,m)-s
end do
! sitten U:n vaakarivi
U(m,m)=1
do j=m+1,n
s=0.0
do k=1,m-1
s=s+L(m,k)*U(k,j)
end do
U(j,m)=(A(j,m)-s)/L(m,m)
end do
end do
Tässä voidaan itse asiassa korvata L ja U A:lla, jolloin hajotelma kirjoitetaan alkuperäisen matriisin päälle.
Hajotelman avulla lausuttuna yhtälöryhmä Ax=b tulee muotoon
LU x = b.
Tämä vastaa yhtälöitä
L y = b,
U x = y.
Kummassakin kerroinmatriisi on kolmiomatriisi, joten yhtälöryhmä on helppo ratkaista pelkällä takaisinsijoituksella.
Aluksi ratkaistaan vektori y yhtälöryhmästä
L y = b.
l11 0 ... 0 y1 b1 l21 l22 ... 0 y2 b2 . . . . . = . . . . ln1 ln2 ... lnn yn b1eli
l11 y1 = b1,
l21 y1+l22 y2 = b2,
.
.
.
ln1 y1+ln2 y2 + ... + lnn yn = bn.
y1 = b1/l11,
y2 = (b2-l21 y1)/l22,
.
.
.
yn = (bn-l_n1 y1-ln2 y2 - ... - ln,n-1yn-1)/lnn.
Lukua yi laskettaessa ei enää
tarvita aikaisempia arvoja bj, j < i,
joten ratkaisu voidaan kirjoittaa
b-vektorin päälle.
Seuraavaksi ratkaistaa yhtälö U x = y:
1 u12 ... u1n x1 y1 0 1 ... u2n x2 y2 . . . . = . . 0 0 ... 1 xn yneli
x1 + u12x2 + ... + u1nx1 = y1, . . . xn-1 + un-1,n xn = yn-1, xn = yn.
xn = yn, xn-1 = yn-1-un-1,n xn, . . . x1 = y1- u12x2 - ... - u1nxn.Taaskin ratkaisut voidaan kirjoittaa vakioiden yi tilalle.
Esimerkki
1 2 2 x1 1 1 1 0 x2 = 2 2 -1 1 x3 3Edellä lasketun LU-hajotelman avulla tämä vastaa yhtälöryhmää
1 0 0 1 2 2 x1 1 1 -1 0 0 1 2 x1 = 2 2 -5 7 0 0 1 x1 0
y1 = 1/1 = 1, y2 = (2-1)/(-1)= -1, y3 = (0-2-5)/7 = -1.
x3 = -1, x2 = -1+2 = 1, x1 = 1-2+2 = 1.
Menetelmä on kätevä, jos ratkaistavana on useita yhtälöryhmiä, joissa esiintyy sama kerroinmatriisi. Matriisin LU-hajotelma tarvitsee laskea vain kerran.
M x = (M - A) x + b,
missä M on mielivaltainen matriisi. Jos nyt xi on kierroksella i laskettu ratkaisun approksimaatio, seuraava iteraatti saadaan kaavasta
M xi+1 = (M - A) xi + b.
Tämä on siis lineaarinen yhtälöryhmä, josta xi+1 on ratkaistava. Jotta menetelmä olisi mielekäs, matriisin M on luonnollisesti oltava sellainen, että tämän yhtälöryhmän ratkaiseminen on oleellisesti helpompaa kuin alkuperäisen yhtälöryhmän.
Yksinkertaisin vaihtoehto on valita M-matriisiksi yksikkömatriisi. Ratkaisuvektorin uusi iteraatti on silloin suoraan oikean puolen vektori.
Lähes yhtä yksinkertainen menetelmä saadaan, jos M-matriisiksi valitaan alkuperäisen kerroinmatriisin A} lävistäjä. Tämä tunnetaan Jacobin menetelmänä.
Menetelmä suppenee kaikille alkuarvoille, jos kerroinmatriisin diagonaali on dominoiva:
| aii | > sum_{j \ne i} | aij |, i=1, ..., n.
program jacobi
! Lineaarisen yhtaloryhman
! ratkaisu Jacobin iteraatiolla
integer, parameter :: n=3
real, parameter :: epsilon = 0.00001
integer i, j, iter
real d, s, x(n), y(n), A(n,n), M(n), b(n)
A(1,:) = (/ 3, 1, 1 /)
A(2,:) = (/ 1, -3, -1 /)
A(3,:) = (/ 2, 1, -4 /)
b = (/1, 2, 2/)
! M-matriisi; vain lavistajaalkiot tarvitaan
do i=1,n
M(i)=A(i,i)
end do
! A korvataan A-M:lla
do i=1,n
A(i,i)=0
end do
! ratkaisuvektorin alkuarvo
x=0.0
! iteroidaan, kunnes perakkaiset approksimaatiot
! x ja y eroavat toisistaan riittavan vahan
d=1.0
iter=0
do while (d > epsilon)
! iteraatiokaavan oikean puolen matriisi
do i=1,n
s=0.0
do j=1,n
s=s-A(i,j)*x(j)
end do
y(i) = (s+b(i))/M(i)
end do
! perakkaisten iteraattien ero
d=0.0
do i=1,n
d=d+(x(i)-y(i))**2
end do
! paivitetaan ratkaisuvektori
x=y
write (6,*) x, d
iter = iter+1
if (iter > 10) exit
end do
end
0.3333333 -0.6666667 -0.5000000 0.8055556
0.7222223 -0.3888889 -0.5000000 0.2283951
0.6296296 -0.2592592 -0.2361111 9.5014609E-02
0.4984568 -0.3780864 -0.2500000 3.1519126E-02
0.5426955 -0.4171811 -0.3452932 1.2566250E-02
0.5874915 -0.3706704 -0.3329476 4.3223356E-03
0.5678726 -0.3598537 -0.2989219 1.6596466E-03
0.5529252 -0.3777352 -0.3060271 5.9365941E-04
0.5612541 -0.3803492 -0.3179712 2.1886613E-04
0.5661068 -0.3735916 -0.3144603 8.1541584E-05
0.5626839 -0.3731443 -0.3103445 2.8855728E-05
Skaalaus. Jakajaksi saattaa tulla hyvin pieniä tai hyvin suuria lukuja, jolloin kertoimien suuruusluokka voi muuttua merkittävästi.
Suppeneminen. Iteratiivisissa menetelmissä on seurattava ratkaisun suppenemista. Varsinaisen suppenemiskriteerin lisäksi mukana on hyvä olla jokin rajoitus kierrosmäärälle.
Yli- ja alivuodot. Laskutoimituksen seurauksena voi syntyä luku, joka on liian suuri esitettäväksi tietokoneessa. Ongelma yleensä vältetään sopivalla skaalauksella.
Esimerkki: Kahanin yhtälöpari:
1.2969 x + 0.8648 y = 0.8642 0.2161 x + 0.1441 y = 0.1440Kun lasketaan 8 desimaalilla, saadaan ratkaisu
x = 1.33317912, y=-1.0
Residuaalin
r = A x - b
komponentit ovat luokkaa 10-8, siis sama kuin laskentatarkkuus.
Oikea ratkaisu on kuitenkin x=2, y=-2!
cond(A) = || A || || A-1 ||.
Yhtälöryhmää ratkaistaessa lähtöarvojen virhe kasvaa suunnilleen kertoimella cond(A). Helposti laskettava matriisinormi on
|| A || _infty = max_i sum_j | aij |.
A = 1.2969 0.8648
0.2161 0.1441
A-1 = 108 x 0.1441 -0.8648
-0.2161 1.2969
|| A || \approx 2 ja
|| A-1 = 2 108, joten
cond}(A) \approx 4 108.
det A = det LU = det L det U.
L ja U ovat kolmiomatriiseja, joiden determinantti on lävistäjällä olevien alkioiden tulo. Matriisin U lävistäjällä on pelkkiä ykkösiä, joten sen determinantti on 1.
det A = l11l22 ... lnn = prod_i=1 .. n lii.
A A-1 = A-1A = I.
Käänteismatriisia ei pidä laskea, ellei sitä todellakin tarvita johonkin. X on haluttu käänteimatriisi, jos
A X = I.
Tämä vastaa yhtälöitä
A xi = ei, i=1, ... ,n,
missä xi on matriisin X i:s pystyrivi ja yksikkömatriisin ei i:s pystyrivi.