KURZUS: Számítási módszerek

MODUL: Matematikai számítások MATLAB-bal

5. lecke: Lineáris algebrai alkalmazások

Cél: Ahogy már az Excel számítások modulban említettük, a lineáris egyenletrendszerek témaköre a gyakorlati mérnöki-gazdasági alkalmazások különösen gazdag és szerteágazó volta miatt - tananyagunkban is - kiemelt fontossággal bír. Mivel a MATLAB rendkívül impozáns és hatékony megoldó-apparátust tud nyújtani ehhez a feladattípushoz, feltétlenül érdemes ezt a támogatást megismernünk és használnunk. Szintén kitüntetett szerephez jut a mérnöki tudományokban a sajátérték-sajátvektor problémakör - és itt is nagyon jó a MATLAB-os támogatás. Ebben a leckében már az eddigi tudás egyfajta gyakorlati szintézisét is megcélozzuk, mintegy "feltéve a koronát" az elméletre (a fáradtságos munkával megtanult mátrixok).

Az Excel modulban már megismert fogalmakra építve ez a lecke viszonylag könnyen elsajátítható, még azzal a megjegyzéssel is, hogy - természetesen - nem tudunk mindent "mechanikusan" átvenni a táblázatkezelős környezetből. A leckében új anyagként szereplő további részek - a tananyagban előírt szinten - szintén elég egyszerűen megtanulhatók.

A lecke elsajátítása után - a probléma megfelelő elemzése és értelmezése után - Ön is képes lesz a saját munkájában a módszerek gyors és ügyes alkalmazására. Az itt megszerzett tudás sok konkrét gyakorlati feladatoknál kamatoztatható!

Követelmények: Ön akkor sajátította el megfelelően a tananyagot, ha (a MATLAB segítségével)

  • elő tudja állítani két vektor skalárszorzatát, vektoriális szorzatát, hajlásszögét;
  • meg tudja határozni egy mátrix és egy vektor egyes, kettes és végtelen normáját;
  • ki tudja számolni egy (kisebb méretű) mátrix sajátértékeit és sajátvektorait;
  • el tudja dönteni, hogy egy négyzetes lineáris egyenletrendszer megoldása egyértelmű-e;
  • képes az inverz mátrixos és a balosztós módszerrel is az egyértelmű eset megoldására;
  • szinguláris (négyzetes) együttható mátrix esetén el tudja dönteni (valamelyik tanult módszerrel), hogy a rendszer összefüggő vagy ellentmondásos-e;
  • saját szavaival meg tudja fogalmazni, hogy milyen kitüntetett x vektorokat ("megoldás") keresünk végtelen sok, illetve 0 megoldás létezése esetén;
  • képes meghatározni a kitüntetett megoldást az összefüggő és a túlhatározott, ill. ellentmondó esetben;
  • mindhárom fenti esetben tudja ellenőrizni, hogy a kapott megoldások valóban eleget tesznek az előírt feltételeknek;
  • kategorizálni tudja az egyes megoldási módszereket hatékonyság szempontjából.

Időszükséglet: A tananyag elsajátításához (a feladatok megoldásával együtt) hozzávetőlegesen 3 órára lesz szüksége.

Kulcsfogalmak

  • Skalárszorzat, vektoriális szorzat, vektorok hajlásszöge;
  • Egyes, kettes és végtelen norma; sajátérték, sajátvektor;
  • Lineáris egyenletrendszer, együtthatómátrix, ennek rangja és determinánsa;
  • Kibővített rendszer és rangja, eltolt együtthatómátrix és determinánsa;
  • Lineáris egyenletrendszer megoldása és kitüntetett megoldása, eltérési hiba.
  • A megoldási módszer hatékonysága.
Speciális szorzatok

Mérnöki problémák megoldása során gyakran szükséges speciális vektorszorzatok kiszámítása. A MATLAB-ba - elvárásainknak megfelelően - beépítettek megfelelő függvényeket ilyen szorzatok meghatározására. A két legegyszerűbb művelet a következő függvényekkel állítható elő:

  • dot(x, y) - euklideszi skalárszorzatot állít elő (szorzatösszeg, nemcsak 3-dimenzióban használható), ahol x és y azonos méretű vektorok.
  • cross(a, b) - vektoriális szorzat, ahol a és b 3-elemű vektorok, az eredmény az a és b vektorok által meghatározott síkra merőleges.

Megjegyezzük, hogy a skalárszorzatot a sum(x, y) képlet is előállítja.

Egyszerű használati példák:

>> a = [3,4,0], b = [-4,3,0], dot(a, b)  % a vektorok merőlegesek
a =
    3    4    0
b =
    -4    3    0
ans =    0
>> a = [1 0 0], b = [0 1 0], c = cross(a, b)
a =
    1    0    0
b =
    0    1    0
c =
    0    0    1

A fenti példa módosításával állítsa elő két párhuzamos vektor skaláris szorzatát.

Határozza meg az előző példában szereplő a, b és c vektorokra a b×a, a×c és c×a vektoriális szorzatokat.

A skaláris szorzattal ellentétben a vektoriális szorzat meghatározása csak háromdimenziós vektorok esetében megengedett.

>> a = [1 0]; b = [0 1]; cross(a, b)
??? Error using ==> cross at 37
A and B must have at least one dimension of length 3.

A két parancs kombinálva is használható: ha x, y és z 3-elemű vektorok, akkor a dot(cross(x,y),z) kifejezés a három vektor vegyesszorzatát jelenti, és a három vektor által kifeszített paralelepipedon térfogatát adja. Ha ez 0-nak adódik, akkor a 3 vektor lineárisan nem független, azaz egy síkban fekszenek.

Magasabb geometriai ismereteket felhasználva tudjuk, hogy a vegyesszorzat a három vektorból képzett mátrix determinánsának kifejtésével, azaz a det([x; y; z]), ill. det([x, y, z]) parancsokkal is meghatározható.

>> a = rand(1,3), b = rand(1,3), c = rand(1,3)  % véletlen vektorok
a =
    0.8147    0.9058    0.1270
b =
    0.9134    0.6324    0.0975
c =
    0.2785    0.5469    0.9575
>> dot(cross(a, b), c)  % vegyesszorzat a definíció szerint
ans =  -0.2767
>> det([a; b; c])  % vegyesszorzat determinánsszámítással
ans =  -0.2767

Próbálja ki, hogy a determinánsos vegyesszorzat-meghatározás ugyanazt az eredményt adja akkor is, ha az a, b és c vektorok transzponáltjait (oszlopvektorok) egymás mellé helyezve építjük fel az alapmátrixot.

Vektorok hajlásszöge a MATLAB-ban a subspace(x, y) paranccsal határozható meg, az eredmény a két oszlopvektor "szöge" radiánban. Sorvektorokra nem kapunk helyes eredményt! (A parancs 3-nál magasabb dimenzióban is használható.)

>> subspace(a', b')  % oszlopvektorok: helyes eredmény (radiánban)
ans =    0.2322
>> subspace(a, b)  % sorvektorok: az eredmény hibás
ans =    0

A koszinusztételből tudjuk, hogy a hajlásszög a skaláris szorzatból kiindulva is meghatározható. Így a subspace parancs helyett használhatjuk az acos(dot(x, y)/(norm(x)*norm(y))) kifejezést, ahol a norm(x) = sqrt(dot(x, x)) kifejezés az euklideszi norma, azaz az x vektor "hossza" (a koordináták négyzetösszegéből vont gyök, lásd lent), és ez képlet már sorvektorokra is működik.

>> acos(dot(a, b)/(norm(a)*norm(b)))  % sorvektorok
ans =    0.2322

Nullvektor esetében természetesen az eredmény értelmetlen lenne, de mint tudjuk, ilyenkor a vektor iránya sem értelmezett, tehát már a kérdésfeltevés is hibás.

>> d = [0 0 0];  % nullvektort adunk meg
>> cross(a, d)  % a vektoriális szorzat nullvektor
ans =
    0    0    0
>> acos(dot(a, d)/(norm(a)*norm(d)))  % a hajlásszög nem értelmezett
ans =  NaN

Kockás füzetébe rajzoljon le két olyan síkvektort, amelyek egymással 45 fokos szöget zárnak be. Ezeket a vektorokat vegye fel a MATLAB-ba x és y azonosítóval. Ellenőrizze, hogy az arkuszkoszinuszos képlet felhasználásával valóban megkapja a 45 fokos bezárt szöget.

Tipp: a fokban megkapandó válaszhoz használja az acosd függvényt!

Megjegyezzük még, hogy a subspace parancs mátrix paraméterekkel is használható. Ilyenkor a subspace(A, B) függvény az A és B mátrixok oszlopvektorai által meghatározott alterek szögét számolja ki. Ennek akkor van értelme, ha az A és B mátrixok oszlopainak száma kisebb, mint a sorok száma, azaz valódi alterekről van szó. További részletek a súgóban olvashatók.

Vektornormák

A vektorok hajlásszögének meghatározásánál találkoztunk a norma fogalmával. A következőkben áttekintjük a meghatározás lehetőségeit és a rendelkezésre álló MATLAB-os támogatást.

Egy vektor normája a vektor hosszának több dimenzióban való általánosítása. A vektornorma tulajdonságai - definíció szerint - a következők:

  • Nem negatív, és csak akkor 0, ha a vektor maga a nullvektor;
  • A vektor megnyújtottjának (konstansszorosának) normáját úgy kapjuk, hogy a vektor normáját megszorozzuk a nyújtási tényező abszolút értékével;
  • Érvényes a háromszög-egyenlőtlenség.

Ezt a definíciót a szám n-esek vektorterében többféle módon is lehet teljesíteni, ezért a matematikában többféle normát definiáltak. Ezek konvergens vektorsorozatokra ekvivalensek, azaz ha egy xk vektorsorozat eltérése egy x0 vektortól valamelyik normában mérve nullához konvergál, akkor a többi normában mérve is.

A norma fogalma analóg módon kiterjeszthető a mátrixokra is (vektorokból képzett vektort tekintve).

A MATLAB segítségével többféle normát is használhatunk, ezek közül mi az egyes, a kettes és a végtelen normáról teszünk említést (így például nem tárgyaljuk a MATLAB által szintén ismert Frobenius-féle normát).

>> help norm  % a válasz részletét közöljük
NORM  Matrix or vector norm.
    For matrices...
      NORM(X) is the 2-norm of X.
      NORM(X,2) is the same as NORM(X).
      NORM(X,1) is the 1-norm of X.  % max(sum(abs(A)))
      NORM(X,inf) is the infinity norm of X.  % max(sum(abs(A')))
...
    For vectors...
      NORM(V,P) = sum(abs(V).^P)^(1/P).
      NORM(V) = norm(V,2).
      NORM(V,inf) = max(abs(V)).
      NORM(V,-inf) = min(abs(V)).

Láthatjuk, hogy mátrixokra az egyes norma a legnagyobb oszlopösszeg, a végtelen norma pedig a legnagyobb sorösszeg. A kettes normát - ez az alapértelmezett norma - mátrixokra úgy kapjuk meg, hogy az X'*X mátrix abszolútértékben legnagyobb sajátértékéből gyököt vonunk (a sajátérték fogalmát lásd a következő alfejezetben):

>> A = [1 3 6; 2 0 -1; 6 1 2]  % tesztmátrix
A =
    1    3    6
    2    0    -1
    6    1    2
>> norm(A, 1)  % legnagyobb oszlopösszeg
ans =    9
>> norm(A, inf)  % legnagyobb sorösszeg
ans =    10
>> norm(A)  % alapértelmezés: kettes norma
ans =    8.0548
>> sqrt(max(abs(eig(A'*A))))  % ugyanezt kapjuk a képlettel is
ans =    8.0548

Hozzon létre egy 3×3-as, véletlen egészeket tartalmazó tesztmátrixot (elemei legyenek 1 és 19 közötti egész számok). Próbálja ki ezen a mátrixon az egyes, a kettes és a végtelen normát!

Vektorokra a kettes norma (ez az alapértelmezés) a hagyományos euklideszi távolságot, adja, a végtelen norma az abszolút értékben legnagyobb koordináta, az egyes norma pedig a "sima" koordinátaösszeg.

>> a = [0 1 2 3 4]
a =
    0    1    2    3    4
>> norm(a)  % kettes norma, koordináták négyzetösszegéből vont gyök
ans =    5.4772
>> norm(a, 1)  % egyes norma, koordinátaösszeg
ans =    10
>> norm(a, inf)  % végtelen norma, most: legnagyobb érték
ans =    4

A normák jobb megismeréséhez természetesen ajánljuk még az egyetemeken folyó lineáris algebra kurzusok lelkiismeretes látogatását, illetve a megfelelő jegyzetek (Kiss-Krebsz 2005) tanulmányozását.

Sajátérték, sajátvektor

A műszaki matematika, statisztika stb. nagyon sok problémájának megoldása az ún. sajátérték-feladat, vagy az általánosított sajátérték-feladat megoldására támaszkodik. Ez a feladat adódik például egy híd terheléspróbájából megnyúlások alapján a híd rezonancia-frekvenciáinak kiszámítására. A sajátérték-feladat legegyszerűbb megfogalmazása:

Au = λu

Itt azt kell tudni, hogy ha az A lineáris transzformációt az A mátrix reprezentálja, azaz A(x) = A*x, akkor melyek azok az irányok (u sajátvektorok), amelyeket az A transzformáció csak megnyújt? Hasonlóképpen milyen nyújtási tényezők (sajátértékek) adódnak? Ilyen feladat lehet például egy merev test esetében annak a forgástengelynek (szabad tengely) megállapítása, amely körüli forgásnál eltűnnek a deviációs momentumok.

Ahogy azt matematikából (lineáris algebra) tudjuk, az A négyzetes mátrix sajátértékeit a karakterisztikus polinomjának zérushelyei adják. A karakterisztikus polinomot az EλA mátrix determinánsának kifejtésével tudjuk felírni, ezután lehetne a zérushelyeket meghatározni.

A fentiek alapján hány sajátértéke lehet egy n×n-es mátrixnak?

Ezt persze nekünk a MATLAB-ban nem kell így kiszámolnunk, mert kényelmes beépített eszközök állnak rendelkezésünkre. A sajátérték feladat a legegyszerűbb módon az eig parancs használatával oldható meg.

Ha nem rendelkezünk másként, akkor a függvény egy oszlopvektorban a megadott mátrix sajátértékeit adja válaszul.

>> A = [3, 4, 1; 2, 0, 2; -1, 1, 1]
A =
    3    4    1
    2    0    2
    -1    1    1
>> eig(A)  % sajátértékek
ans =
    4.3166
  -2.3166
    2.0000

A használat másik lehetséges esete, hogy a visszatérési értékek számára paraméterként egy mátrixpár-tárolót adunk meg (például: U és Lambda). Ekkor két új n×n-es mátrix keletkezik, az első az 1-re normált sajátvektorokat fogja tartalmazni (U), a második pedig - főátlójában - a megfelelő sajátértékeket (Lambda).

>> [U Lambda]= eig(A)  % a sajátvektorokat is meghatározunk
U =
    0.9218    0.5074  -0.7071
    0.3468  -0.7707  -0.0000
  -0.1734    0.3854    0.7071
Lambda =
    4.3166        0        0
        0  -2.3166        0
        0        0    2.0000

Ellenőrzésként megnézzük, hogy az A mátrix valóban csak nyújtja az U-beli vektorokat:

>> A*U
ans =
    3.9790  -1.1754  -1.4142
    1.4968    1.7855  -0.0000
  -0.7484  -0.8928    1.4142

Jól látható, hogy ugyanazt kapjuk, mintha a sajátvektorokat rendre csak a megfelelő (!) sajátértékekkel szoroznánk össze:

>> U*Lambda
ans =
    3.9790  -1.1754  -1.4142
    1.4968    1.7855  -0.0000
  -0.7484  -0.8928    1.4142

A fenti, véletlen egészeket tartalmazó 3×3-as tesztmátrixszal próbálja ki a sajátértékek és sajátvektorok meghatározását! Ellenőrizze, hogy az eredeti mátrix valóban csak nyújtja-e a kapott sajátvektorokat!

Megjegyezzük, hogy a sajátérték probléma általánosabb vizsgálatához és megoldásához (nagy mátrixok, ritka mátrixok stb.) további, mélyebb ismeretekre van szükség. Ezt mi a Számítási módszerek tananyagban terjedelmi okokból nem tudjuk tárgyalni. Az érdeklődők minderről további részleteket olvashatnak Stoyan Gisbert MATLAB című könyvében.

Kondíciószám

Egy X négyzetes és invertálható mátrix kondíciószáma azt mutatja meg, hogy egy mátrixalgebrai feladat megoldása mennyire érzékeny az X mátrix pici megváltoztatására. Ha a mérnöki számítások során figyelmen kívül hagyjuk azt, hogy egy mátrix rosszul kondicionált, az nagyon kellemetlen következményekkel járhat; ezért a jelenséget fontos ismerni (és kezelni tudni).

A definíció a következő:

cond(X,p)= X p X 1 p aholp=1,2,inf,'fro'

Ez a definíció a 2-es normában a legnagyobb és legkisebb szinguláris érték hányadosa (X mátrix szinguláris értékeinek az X'*X mátrix - ami egyébként már szimmetrikus mátrix lesz - sajátértékeit nevezzük, lásd még help svd).

>> cond(A)  % kondíciószám (kettes normában)
ans =    3.3188
>> norm(A)*norm(inv(A))  % a definíció szerint felírva
ans =    3.3188
>> sqrt(max(abs(eig(A'*A))))/sqrt(min(abs(eig(A'*A)))) % és mint a...
ans =    3.3188
>> max(svd(A))/min(svd(A))  % max és min szing. értékek hányadosa
ans =    3.3188

A nagy kondíciószám közel szinguláris mátrixot jelent. A MATLAB-ban az rcond(X) reciprok kondíciószámot is lehet használni, ami az 1-es normával számolt kondíciószám reciproka. Jól kondicionált mátrixokra az rcond értéke 1-hez közeli, míg egy rosszul kondicionált mátrix esetén az rcond értéke eps nagyságrendű.

>> rcond(A)  % beépített függvény
ans =    0.1667
>> 1/(norm(A,1)*norm(inv(A),1))  % mi kiszámoljuk ki
ans =    0.1667

A MATLAB modul 3. leckéjében már említettük - és szemléletesen meg is mutattuk -, hogy a Hilbert-mátrixok gyengén kondicionált mátrixok. Így rájuk a cond és az rcond parancs nagy, illetve igen kicsi értéket szolgáltat.

>> cond(hilb(10)), rcond(hilb(10))
ans =    1.602466539329798e+013
ans =    2.828556315915063e-014

A fent bemutatott példák alapján ellenőrizze, hogy ugyanezeket a válaszokat kapjuk akkor is, ha a cond és az rcond értékeket a definíció alapján számoljuk ki.

Lineáris egyenletrendszerek megoldása - egyértelmű eset

Idézze fel az Excel számítások modul 2. leckéjében tanultakat a lineáris egyenletrendszerekkel kapcsolatos fontos fogalmakról és a megoldás lehetőségeiről.

Már korábban láttuk (Excel számítások modul, 2. lecke), hogy egy - négyzetes - Ax = b lineáris egyenletrendszer A együtthatómátrixának és az inhomogén jobboldali b oszlopvektor ismeretében az x megoldásvektor egyszerűen számítható. A megoldhatóság feltétele az, hogy az A mátrix ne legyen szinguláris, azaz a determinánsa ne legyen 0. Mint tudjuk a MATLAB modul 2. leckéjéből, ugyanez a feltétel a rang vizsgálatával is mérhető.

>> A = [3 4 1; -1 2 5; 7 -2 3]
A =
    3    4    1
    -1    2    5
    7    -2    3
>> b = [18; 8; 32]
b =
    18
    8
    32
>> det(A)  % nem nulla determináns, egyértelmű megoldás létezik
ans =  188
>> rank(A), rank(A) == length(A)  % a rang is ezt mutatja
ans =    3
ans =    1

Megoldjuk az egyenletrendszert az inverz mátrixos módszerrel (inv parancs; az inv helyett itt a pinv parancs is használható, részletesen lásd a következő alfejezetben), majd ellenőrizzük a megoldást:

>> x = inv(A)*b
x =
    4.0000
    1.0000
    2.0000
>> A*x  % ellenőrzés 1. - valóban b-t kapjuk
ans =
    18
    8
    32
>> A*x == b  % ellenőrzés 2. - a koordináták mind megegyeznek
ans =
    1
    1
    1

A MATLAB modul 2. leckéjében azt is tanultuk, hogy a MATLAB-ban a balosztó művelet segítségével A-hoz és b-hez A\b-vel előállítható egy olyan x mátrix, amelyre Ax = b. Ezt a módszert is felhasználhatjuk a lineáris egyenletrendszer egyértelmű megoldásához:

>> x = A\b
x =
    4.0000
    1.0000
    2.0000

Látható, hogy az egyértelmű eset megoldásához mindkét módszer egyaránt használható. Arra a természetesen felmerülő kérdésre, hogy a kettő közül melyik módszer a hatékonyabb (gyorsabb), később térünk ki.

Megjegyezzük, hogy kis numerikus hibák miatt előfordulhat, hogy az A*x visszaszorzásnál nem kapjuk meg pontosan a b vektort (csak kisebb hibával), ill. az A*x == b ellenőrzés nem ad logikai igaz választ (egyes koordinátákra, vagy egyre sem). Ez azonban természetesen nem jelenti azt, hogy hibásan dolgoztunk; a jelenséget a numerikus rendszerek természetes "velejárójaként" kell elkönyvelnünk.

A fenti példa alapján ellenőrizze, hogy az A = [4 5 6; 5 6 8; 1 7 7], b = [23; 30; 11] lineáris egyenletrendszer megoldása egyértelmű. Oldja meg a rendszert az inverz mátrixos és a balosztós módszerrel! Ellenőrizze a megoldást a tanult módon! Mit ad az egyes esetekben az A*x == b ellenőrzés? Hogyan magyarázza a jelenséget?

Lineáris egyenletrendszerek - ellentmondásos eset

Idézze fel az Excel számítások modul 2. leckéjében tanultakat az összefüggő és az ellentmondásos rendszerek megoldásáról, illetve a MATLAB modul 2. leckéjében tanultakat a pszeudoinverz meghatározásáról!

Módosítsuk az előző A mátrixot a következő módon:

>> A(3, 2:3) = [6 -3]
A =
    3    4    1
    -1    2    5
    7    6    -3

Ha most az inverz mátrixos módszert választjuk a megoldáshoz (pontosabban: a megoldási kísérlethez), az eredmény megjelenítésénél a következő üzenetet kapjuk:

>> x = inv(A)*b  % x = A\b is ugyanezt adná!
Warning: Matrix is close to singular or badly scaled.
        Results may be inaccurate. RCOND = 8.410780e-018.
x =
  1.0e+015 *
    8.1065
  -7.2058
    4.5036

Látjuk, hogy a MATLAB figyelmeztet bennünket arra, hogy a mátrix közel szinguláris (gyengén kondicionált), és az rcond értéke igen kicsi, emiatt az eredmény pontatlan lehet (és valóban az is).

Mennyi lesz a példában a cond érték? (A pontos meghatározás előtt végezzen becslést!)

Visszaszorzással ellenőrizze, hogy a kapott x vektor ténylegesen nem megoldása a lineáris egyenletrendszernek.

A problémát az okozta, hogy a módosítás után elfelejtettük ellenőrizni az A mátrix determinánsát, illetve rangját (ez a lépés pedig nem hagyható ki). Némi elemzés után látható, hogy ez esetben a mátrix harmadik sora az első sor kétszeresének és a második sornak a különbsége lesz, így a sorvektorok (és persze az oszlopvektorok is) összefüggő rendszert alkotnak. A mátrix tehát valóban szinguláris és így ténylegesen nem is invertálható (tehát pusztán a "szerencsének" köszönhettük, hogy egyáltalán kaptunk Inf értékektől különböző választ):

>> det(A), rank(A)
ans =    0
ans =    2

Az Excel számítások modul 2. leckéjében már említettük, de itt is újra megjegyezzük azt, hogy előfordulhat - kis kerekítési pontatlanságok miatt - a következő jelenség: a determináns ténylegesen 0 ugyan, de az ellenőrzés nem pontosan nullát, hanem egy nullához nagyon közeli eredményt ad (ekkor pedig az értelmezéshez még tovább kell gondolkodnunk). Ezért jobb eszköz a vizsgálathoz a rang: ilyenkor azonnal biztos választ kapunk.

Ahogy az Excel számítások modul 2. leckéjében már tanultuk, a továbbiakban el kell döntenünk, hogy a teljes rendszer összefüggő vagy ellentmondásos-e. Az Excelben ilyenkor az eltolt rendszer determinánsának vizsgálatával tudtunk döntést hozni.

Természetesen a MATLAB-ban is dolgozhatunk így; azonban a MATLAB a rang vizsgálatával ilyen téren "többet tud", mint az Excel, ezért más módszerek is rendelkezésre állnak.

A legkézenfekvőbb a kibővített rendszer rangjának ellenőrzése:

>> rank([A b])  % a kibővített rendszer rangja
ans =    3

A kibővített rendszer vizsgálatával tehát arra a következtetésre jutunk, hogy a teljes rendszer ellentmondásos.

Az eltolt rendszer vizsgálatával is ugyanígy célhoz érünk:

>> eltolt = [A(:, 2:3) b]  % eltolt rendszer felépítése
eltolt =
    4    1    18
    2    5    8
    6    -3    32
>> det(eltolt), rank(eltolt)  % vizsgálat
ans =    72
ans =    3

Ahogy már szintén tanultuk (Excel számítások modul, 2. lecke), az ellentmondásos esetben pontos megoldás nem állítható elő, de van lehetőség egy valamilyen szempontból optimálisnak tekinthető, közelítő megoldás meghatározására (az A*x - b hibavektor normája legyen minimális; ez a Moore-Penrose-féle kváziinverz fogalmához vezet).

Mint már tudjuk, az Excellel ellentétben a kváziinverz kiszámítása a MATLAB rendszerben direkt módon támogatott (pinv parancs).

>> x = pinv(A)*b  % az ellentmondásos eset "megoldása"
x =
    2.2392
    2.9725
    0.7255
>> A*x  % ellenőrzés
ans =
  19.3333
    7.3333
  31.3333

Meghatározzuk az eltérést:

>> elteres = A*x - b
elteres =
    1.3333
  -0.6667
  -0.6667

Látjuk, hogy a MATLAB a hibák nagyjából egyenletes szétosztásával közelítette az elméletileg nem létező, "elképzelt ideális" megoldást.

Tekintsük a következő egyenletrendszert: A = [1 2; 3 c] és b = [3; 8]. Válassza meg a c paraméter értékét úgy, hogy az A mátrix rangja 1 legyen! (Útmutató: a második sor legyen az első háromszorosa.)

Ellenőrizze, hogy az így előállt rendszer ellentmondásos! Oldja meg a rendszert a tanult módon, majd visszaszorzás után határozza meg a hibavektort! Figyelje meg, hogy a MATLAB nagyjából egyenletesen osztotta el a hibát a koordináták között.

Megjegyezzük, hogy az ellentmondásos eset itt bemutatott típusában (n egyenlet, az A mátrix rangja n - 1, a kibővített rendszer rangja viszont n) a balosztós módszer általában nem használható a megoldáshoz (ahogy fent is láttuk). (A túlhatározott rendszernél viszont alkalmazható lesz; lásd lent.)

Lineáris egyenletrendszerek - összefüggő eset

Módosítsuk most a b vektort is a következő módon:

>> b(3) = 28
b =
    18
    8
    28

Ezzel a teljes rendszer összefüggővé vált, hiszen most már a b vektorra is érvényes, hogy a harmadik sor az első sor kétszeresének és a második sornak a különbsége. Az A mátrix továbbra is szinguláris lesz, viszont így a kibővített rendszer rangja is 2 marad (illetve az eltolt rendszer rangja 2, determinánsa 0):

>> rank([A b])
ans =    2
>> eltolt = [A(:, 2:3) b]
eltolt =
    4    1    18
    2    5    8
    6    -3    28
>> det(eltolt), rank(eltolt)
ans =    0
ans =    2

A vizsgálattal tehát arra a következtetésre jutunk, hogy a rendszer összefüggő (kevesebb független egyenletünk van, mint ismeretlen); ilyenkor végtelen sok megoldás van. Ahogy az Excel számítások modul 2. leckéjében már tanultuk, ekkor a végtelen sok megoldás közül egy speciálisat meghatározunk (kitüntetett megoldás). Azt a megoldásvektort szokták kitüntetettnek választani, amelyiknek a normája minimális (azaz: a "legrövidebb" vektorhosszú; koordinátáinak négyzetösszege minimális).

Megjegyezzük, hogy mi - az Excel számítások modullal megegyezően - az Informatika II. tárgy keretein belül csak azzal az esettel foglalkozunk, amikor a független egyenletek száma csak eggyel kevesebb, mint az ismeretlenek n száma (tehát a rendszer rangja n - 1).

A megoldás a MATLAB-ban most is a pszeudoinverz segítségével állítható elő.

>> x = pinv(A)*b
x =
    1.9882
    2.7882
    0.8824
>> A*x
ans =
  18.0000
    8.0000
  28.0000

Kiszámoljuk a kapott megoldás normáját:

>> norm(x)
ans =    3.5364

Az ellentmondásos esettől eltérően most a balosztós módszerrel is kapunk korrekt megoldást (igaz, a rendszer rögtön egy feltűnő figyelmeztetéssel indít):

>> x = A\b
Warning: Matrix is close to singular or badly scaled.
        Results may be inaccurate. RCOND = 8.410780e-018.
x =
    7.6000
  -2.2000
    4.0000
>> A*x  % visszaszorzással ellenőrizzük, hogy valóban megoldást kaptunk
ans =
  18.0000
    8.0000
  28.0000

Ennek a megoldásnak azonban a normája nagyobb, mint az előzőnek, tehát a megadott feltételek szerint nem optimális:

>> norm(x)
ans =    8.8657

A fenti A = [1 2; 3 6] és b = [3; 8] egyenletrendszerben módosítsa a b vektor második koordinátáját úgy, hogy a rendszer teljesen összefüggő legyen! (Azaz rank([A b]) == 1.)

Oldja meg a rendszert a pszeudoinverzes és a balosztós módszerrel is! Ellenőrizze, hogy valóban mindkét esetben megoldást kapott! Hasonlítsa össze a kapott megoldások normáit!

Megjegyezzük, hogy ugyanez a megoldási módszer alkalmazható abban az esetben is, ha az egyenletrendszer (n - 1) × n-es: n változót és n - 1 független egyenletet tartalmaz (lásd: leckezáró feladatok).

Lineáris egyenletrendszerek - túlhatározott eset

Tekintsük most a következő egyenletrendszert:

>> load At.dat, load bt.dat, At
At =
    1    5    13
    2    3    12
    -4    7    3
    3    1    4
    -2    4    1
>> bt'  % oszlopvektor helyett sorvektorként írattuk ki
ans =
    6    5    3    -1    2

Ez az egyenletrendszer túlhatározott (több az egyenlet, mint az ismeretlen). Ahogy az Excel számítások modul 2. leckéjében már említettük, ilyen rendszer például akkor adódhat, ha a mérési adatokon alapuló egyenletek közé pontatlan mérési eredmények kerülnek (és a megoldásnál szintén a hibakiegyenlítés módszerei jöhetnek szóba).

A túlhatározott eset logikailag az ellentmondó kategóriába sorolható, de persze a fent bemutatott példától eltér a szituáció annyiban, hogy az A mátrix nem négyzetes, így determináns nem számolható, a rang pedig legfeljebb a kisebb méret lehet.

>> rank(At)
ans =    3

A kibővített rendszer vizsgálata mutatja, hogy a rendszerben valóban ellentmondás van:

>> rank([At bt])
ans =    4

A kitüntetett megoldás (amikor az A*x - b hibavektor normája minimális) most is előállítható a pszeudoinverzzel.

>> pinv(At)  % a pszeudoinverz kiíratása
ans =
  -0.0618  -0.0690  -0.0331    0.4118    0.0841
  -0.0434  -0.0807    0.0932    0.2802    0.1315
    0.0614    0.0735  -0.0235  -0.1382  -0.0562
>> xt = pinv(At)*bt
xt =
  -1.0591
  -0.4012
    0.6909

Most azonban a balosztó is ugyanezt a megoldást adja:

>> xt = At\bt
xt =
  -1.0591
  -0.4012
    0.6909

Tudjuk, hogy a megoldás nem lehet pontos. Ilyenkor mindig érdemes az eltérésvektort is meghatározni.

>> elteres = At*xt - bt
elteres =
  -0.0830
  -0.0306
    0.5007
    0.1852
  -0.7957

Ahogy már korábban tanultuk, a pszeudoinverz előállítható az eredeti mátrix és transzponáltjának speciális szorzataként is. Konkrétan, ha egy At mátrixnak több sora van, mint oszlopa, és rangja az oszlopok száma, akkor a pszeudoinverz az (At'*At)\At' módon számolódik.

>> (At'*At)\At'  % a fent pinv-vel számolt eredményt kaptuk
ans =
  -0.0618  -0.0690  -0.0331    0.4118    0.0841
  -0.0434  -0.0807    0.0932    0.2802    0.1315
    0.0614    0.0735  -0.0235  -0.1382  -0.0562

Ellenőrizze, hogy a pszeudoinverz ekkor At\eye(sorok száma) módon is számolható!

Megjegyezzük még, hogy mi a Számítási módszerek tárgy keretein belül a túlhatározott rendszereknél is csak olyan esetekkel foglalkozunk, amikor az egyenletek száma csak néhánnyal több, mint az ismeretlenek száma.

Hatékonyságvizsgálat, összefoglalás

Ahogy az egyértelmű eset megoldásánál jeleztük, most visszatérünk arra a problémára, hogy amennyiben több módszer is használható egy lineáris egyenletrendszer megoldására, akkor vajon melyik a hatékonyabb.

Az összehasonlításnál az egyenletrendszer együtthatóit egy speciális tesztmátrix segítségével állítjuk elő (n×n-es Töplitz-mátrix). Itt a célunk az, hogy véletlenszerűnek tekinthető együtthatókat kapjunk, ehhez most 1 és n közötti pozitív egész számokat használunk, megfelelő struktúrába rendezve.

Egy kisebb Töplitz-mátrix a következő:

>> W = toeplitz(5:-1:1)
W =
    5    4    3    2    1
    4    5    4    3    2
    3    4    5    4    3
    2    3    4    5    4
    1    2    3    4    5

Olvassa el a súgóban a toeplitz parancs használatát bemutató részt!

A balosztós és az inverz mátrixos megoldás hatékonyságát fogjuk összehasonlítani az egyértelmű megoldással rendelkező esetben. A T*x = b egyenletrendszert oldjuk meg, ahol T egy 3000×3000-es Töplitz-mátrix, b pedig egy közönséges oszlopvektor. Az időmérést az első MATLAB-os leckében megismert módon (tic-toc páros) oldjuk meg.

>> n = 3000; T = toeplitz(n:-1:1);
>> b = linspace(1,2,n)';  % a mátrixot és a vektort nem íratjuk ki!

A rang ellenőrzésével meggyőződhetünk róla, hogy a rendszer megoldása egyértelmű (a válasz előállítása időigényes, ehelyett érvelhetünk úgy, hogy a konstrukcióból következik a függetlenség, illetve szerezhetünk tapasztalatokat kisebb mátrixokkal).

>> rank(T)  % megkérdezhetjük, de időigényes...
ans =        3000
>> tic; x = T\b; toc
Elapsed time is 0.735710 seconds.
>> tic; x = inv(T)*b; toc
Elapsed time is 9.694462 seconds.
* A végrehajtás egy kétmagos, 2,5 GHz-es gépen történt - gyorsabb négymagos gépeken értelemszerűen rövidebb végrehajtási idő várható; gyengébb gépeken ugyanakkor jóval hosszabb futási is idők is lehetségesek. Ha a példát Ön is ki akarja próbálni így, akkor megfontoltan járjon el!

Láthatjuk, hogy a balosztós módszer nagyjából egy nagyságrenddel gyorsabb.* Gyakorlati példáknál tehát (nagyobb méretű egyenletrendszerek) - az egyértelmű megoldással rendelkező esetekben - feltétlenül érdemes az A\b módszerrel dolgozni.

Kíváncsiak vagyunk még a pinv paranccsal történő megoldás futási idejére is. Itt azonban vigyáznunk kell: nagyméretű, ismeretlen együtthatómátrix esetén ezt a parancsot csak átgondoltan adjuk ki, mert a végrehajtás jóval tovább tart(hat) még az inverz mátrixos esetnél is!

*Ennek oka az, hogy az inverz mátrixos módszernél az LU-felbontást használja fel a MATLAB (ezt mi is néztük a 3. MATLAB-os lecke végén), a pinv-nél pedig egy szinguláris felbontáson alapuló módszert. Mindegyik módszer végrehajtási ideje n 3 -nal arányos, de az utóbbinál a konstans szorzó kb. 35-ször nagyobb! (Gisbert 2005)

Az elmondottak szerint a pinv-es becsült végrehajtási idő a mi gépünkön most 10 perc körül várható, ezért "bevállaljuk".

>> tic; x = pinv(T)*b; toc
Elapsed time is 339.125561 seconds.

A következő táblázatban bemutatjuk a végrehajtási időket a három módszer esetén különböző méretű Töplitz-mátrixokkal.

Nagyméretű lineáris egyenletrendszerek megoldási ideje - összehasonlítás
Méret/Végrehajtási idő (mp) A\binv(A)*bpinv(A)*b
500×5000,0147020,0446241,214775
1000×10000,0664290,37223910,342868
1500×1500 0,1349751,18521134,722225
2000×20000,4495404,34721496,005524
3000×30000,7357109,694462339,125561

1. táblázat

Jól megfigyelhető, hogy az egyes módszerek között a 3000×3000-es mátrixra már megtapasztalt időbeli különbségek itt is jelentkeznek; sőt, alaposabb elemzés finomabb változásokat is kimutathat.

Generáljon egy véletlen egész számokkal feltöltött 500×500-as A együtthatómátrixot a korábban megismert módon (a mátrix elemei 1 és 500 közötti egészek legyenek). A b oszlopvektor legyen pozitív egész számok növekedő sorozata. Oldja meg az így kapott lineáris egyenletrendszert a balosztós, az inverz mátrixos és a pszeudoinverzes módszerrel! Hasonlítsa össze a futási időket!

Ha a gép számolási sebessége lehetővé teszi, növelje a méreteket, és dolgozzon 1000×1000-es, 1500×1500-as stb. együtthatómátrixokkal is.

A következő táblázatban bemutatjuk az alkalmazható megoldási módszereket a tanult esetekben. Megjegyzésként utalunk a módszer gyorsaságára is.

Lineáris egyenletrendszerek megoldási módszerei - összehasonlítás
Módszer/Használható-e inv(A)*bpinv(A)*bA\b
Egyértelmű esetigenigen (de határozottan ez a leglassabb)igen (ez a leggyorsabb)
Végtelen sok megoldás (négyzetes mátrix, ill. (n - 1) × n-es eset)nemigenigen (de a megoldás nem optimális)
Ellentmondásos eset (négyzetes mátrix)nemigennem
Túlhatározott esetnemigenigen (ez a leggyorsabb)

2. táblázat

Összefoglalva, akármelyik esetről is van szó, az Ax = b egyenletrendszer "megoldása" az x = pinv(A)*b módon mindig előállítható, legfeljebb időigényesebb lesz a megoldás.

Önellenőrző kérdések

1. Egy lineáris egyenletrendszer 5×5-ös A együtthatómátrixának determinánsát kiszámítva a következő eredményt kaptuk:

determ =  2.2204e-015
Hány megoldása lehet ez alapján az egyenletrendszernek? Hogyan dolgozna tovább, hogy eldöntse, hogy melyik eset áll fent ténylegesen?

2. Olvassa be MATLAB-ba az A.dat és a b.dat fájlokat!

Az A·x=b lineáris egyenletrendszer együtthatói az A.dat és b.dat fájlban találhatók.

Hány megoldása van az egyenletrendszernek?

Megoldások száma:

Határozza meg az egyenletrendszer minimális normájú x megoldását a Moore-Penrose kváziinverz felhasználásával!

Mennyi az x megoldásvektor-elemek abszolút értékeinek összege?

Az összeg:

Mennyi az x megoldásvektor Euklideszi-normája (elemei négyzetösszegéből vont gyök)?

Az Euklideszi-norma:

Mennyi darab komplex sajátérték van?

A sajátértékek száma:

Mennyi a sajátértékek abszolútértékeinek összege?

Az összeg:

3. A tanult módszerek valamelyikével ellenőrizze, hogy az alábbi táblázattal adott lineáris egyenletrendszer összefüggő!

36-6-7-6 87
-5-754-10 21
-28-41-7 63
4-6-980 6
01-146-23 177

Határozza meg a kitüntetett megoldást! Határozzon meg a MATLAB segítségével egy másik megoldást is! Ellenőrizze, hogy mindkét kapott megoldás valóban megoldás!

4. A tanult módszerek valamelyikével ellenőrizze, hogy az alábbi táblázattal adott lineáris egyenletrendszer ellentmondásos!

36-6-7-6 87
-5-754-10 21
-28-41-7 63
4-6-980 6
01-146-23 176

Határozza meg a kitüntetett "megoldást" (minimalizált eltéréssel)!

5. A tanult módszerek valamelyikével ellenőrizze, hogy az alábbi táblázattal adott 4×5-ös együttható-mátrixú lineáris egyenletrendszer összefüggő!

151314 19
2312-2-3 4
-473-51 -35
31402 14

Határozza meg a kitüntetett megoldást! Határozzon meg a MATLAB segítségével egy másik megoldást is! Ellenőrizze, hogy mindkét kapott megoldás valóban megoldás!

*Tudjuk, hogy amennyiben az As mátrixnak kevesebb sora van, mint oszlopa, és rangja a sorok száma, akkor a pinv() pszeudoinverz az As'/(As*As') módon számolódik. Ellenőrizze ezt a szabályt! Ellenőrizze azt is, hogy a pszeudoinverz ekkor eye(oszlopok száma)/As módon is számítható!