Harmadik labor: Lineáris egyenletrendszerek
Contents
A kondíciószám
Egy mátrix kondíciószámát a cond() függvénnyel kérdezhetjük le, például nem meglepő módon
cond(eye(5),1)
ans =
1
illetve
cond(eye(5),2)
ans =
1
A Matlabban sok nevezetes mátrixcsalád megtalálható beépített függvényként, például a Hilbert-mátrixok is. Ezeknek (i, j)-edik eleme 1/(i+j+1).
H=hilb(4)
H =
1.0000 0.5000 0.3333 0.2500
0.5000 0.3333 0.2500 0.2000
0.3333 0.2500 0.2000 0.1667
0.2500 0.2000 0.1667 0.1429
A Hilbert-mátrixok nevezetesek a magas kondiciószámukról, például
cond(H)
ans = 1.5514e+04
Ugyanakkor a 0 közeli determináns is mutatja, hogy ezek közel szinguláris mátrixok:
det(H)
ans = 1.6534e-07
A magas kondiciószám nem feltétlen jelent 0 közeli determinánst. Vegyük például azt az A mátrixot, amit úgy kapunk, hogy az egységmátrix bal felső elemét kicseréljük egy nagy, pozitív 'a' konstansra. Ekkor det(A)=a és cond(A)=a (hiszen A szimmetrikus, a legnagyobb séje a, a legkisebb 1), például:
A=eye(10); A(1,1)=100; cond(A)
ans = 100
det(A)
ans = 100
Másfelől léteznek 0 közeli determinánssal alacsony kondiciószámú mátrixok. Megkeresésükhöz csak annyit kell kihasználni, hogy a determináns nem skálainvariáns, a kondiciószám viszont igen.
A=eye(10)*0.01; cond(A)
ans =
1
det(A)
ans = 1.0000e-20
Sajátértékek
Az eig() polimorf függvény szolgál a sajátértékek meghatározására:
C=[1 0 3; 4 3 -1; 3 -5 8] eig(C)
C =
1 0 3
4 3 -1
3 -5 8
ans =
-1.5938
4.9175
8.6763
viszont a sajátvektorok lekérdezhetőek:
[Vc,Lc]=eig(C)
Vc =
0.5904 -0.4623 0.3625
-0.6252 -0.6495 0.0920
-0.5105 -0.6037 0.9274
Lc =
-1.5938 0 0
0 4.9175 0
0 0 8.6763
ennek segítségével a spektálfelbontása:
Vc*Lc*inv(Vc)
ans =
1.0000 0.0000 3.0000
4.0000 3.0000 -1.0000
3.0000 -5.0000 8.0000
Mivel C nem szimmetrikus mátrix, ezért V nem ortogonális, így
Vc*Lc*Vc'
ans =
1.6352 2.3543 4.7692
2.3543 1.5251 2.1600
4.7692 2.1600 8.8397
nem a C-t adja vissza. Viszont tudjuk, hogy valós szimmetrikus mátrixokra létezik ortogonális bázis, így például
D=[1 2 3; 2 4 5; 3 5 7]
D =
1 2 3
2 4 5
3 5 7
[V,L]=eig(D)
V =
-0.9008 0.3004 0.3135
-0.0785 -0.8229 0.5628
0.4270 0.4824 0.7648
L =
-0.2478 0 0
0 0.3388 0
0 0 11.9090
V*L*V'
ans =
1.0000 2.0000 3.0000
2.0000 4.0000 5.0000
3.0000 5.0000 7.0000
mert a Matlab normálja V-t, így
V'*V
ans =
1.0000 -0.0000 -0.0000
-0.0000 1.0000 -0.0000
-0.0000 -0.0000 1.0000
de
Vc'*Vc
ans =
1.0000 0.4413 -0.3170
0.4413 1.0000 -0.7872
-0.3170 -0.7872 1.0000
ami mutatja, hogy bár a sajátvektorok egy hosszúak, de nem merőlegesek.
LER-ek megoldása
A linsolve() parancs segítségével oldhatjuk meg a lineáris egyenletrendszereket. Ha nincs speciális tulajdonsága a mátrixnak, akkor a linsolve() Gauss-eliminációt csinál részleges főelemkiválasztással:
A=[1 1 1; 2 4 4; -1 5 8]
A =
1 1 1
2 4 4
-1 5 8
b=[3 ;4 ; 3]
b =
3
4
3
x=linsolve(A,b)
x =
4.0000
-5.0000
4.0000
Az lu() parancs segítségével meg is kaphatjuk a módosított Gauss-elimináció folyamán keletkező P, L, U mátrixokat. Figyeljük meg, hogy L normált mátrix, azaz minden eleme 1-nél kisebb abszolútértékű. A P mátrixból kiderül, hogy volt sorcsere a Gauss-elmináció folyamán.
[L,U,P]=lu(A)
L =
1.0000 0 0
-0.5000 1.0000 0
0.5000 -0.1429 1.0000
U =
2.0000 4.0000 4.0000
0 7.0000 10.0000
0 0 0.4286
P =
0 1 0
0 0 1
1 0 0
Órai feladat
Írjunk olyan függvényt Matlabban, amelynek bemente egy tetszőleges négyzetes mátrix, és végrehajtja a sorcsere nélküli Gauss-eliminációt A-n, a kibővített együtthatómátrixszal most nem foglalkozunk. A kimenet legyen az így kapott L és U mátrix, vagy egy hibaüzenet ha az eljárás elakadt.
function [L, U]=simaGauss(A) U=A; L=eye(size(A,1)); for i=1:size(A,2)-1 if A(i,i)==0 error('A főátlóba 0 került!'); end gausstr=eye(size(A,1)); gausstr(i+1:end,i)= -U(i+1:end,i)./U(i,i); U=gausstr*U; L(i+1:end,i)=-gausstr(i+1:end,i); end