Numerikus AnalízisLineáris egyenletrendszerek megoldásaEgy lineáris egyenletrendszer:
n darab ismeretlen x1,x2, ... ,xn és m darab egyenlet. Az aij együtthatók i=1,2, ... ,m és j=1,2, ... ,n ismert számok és az egyenletek jobb oldalán álló bi, i=1,2, ... ,m számok is azok. Ha n=m, akkor pont annyi ismeretlenünk van, mint egyenletünk. Ekkor van rá reményünk, hogy meg tudjuk oldani az egyenletrendszert. Előfordulhat, hogy az m egyenletből néhány előállítható a többi lineár kombinációjaként, vagy az együtthatók néhány oszlopa állítható elő a többi lineáris kombinációjaként. (Négyzetes mátrixoknál, azaz, ha n=m a két állítás egyenértékű.) Ebben az esetben azt mondjuk, hogy az egyenletrendszer szinguláris, és ekkor nem létezik egyértelmű megoldása. Ha numerikusan próbáljuk megoldani az egyenletrendszert, akkor legalább két problémával kell szembenéznünk:
Gauss--Jordan eliminációEgy mátrix invertálásához a Gauss--Jordan elimináció legalább olyan jó, mint bármely más módszer. Ha egyenletrendszert kívánunk megoldani, akkor a Gauss--Jordan elimináció nem csak az ismeretlenek megoldás vektorát adja eredményül, hanem az együttható mátrix inverzét is. A gyengesége a módszernek,
Az i problémára adhatnánk azt a választ, hogy ha valamilyen új konstans vektorral kívánjuk megoldani az egyenletrendszert, akkor azt csak egyszerűen meg kell szorozni az együtthatómátrix amúgy is adódó inverzével. Ez valóban megtehető, azonban a kerekítési hibák halmozódása miatt ez a megoldás már messze nem lesz olyan pontos, mint azok amelyeket az elimináció során kaptunk meg. Irjunk fel egy négy ismeretlenes egyenletrendszert:
Az ismeretlenek xi, i=1,2,3,4. A konstans tagok bi, i=1,2,3,4. Az együtthatók aij, i=1,2,3,4, j=1,2,3,4. Ugyanaz a megoldás fogja kielégíteni ezt az egyenletrendszert, mint azt, amelyiket ebből úgy kapunk, hogy
Használjuk az egyenletrendszer megoldásához ezt a két műveletet. Osszuk el az első egyenletet a a11-el! Ekkor (1) helyett kapjuk:
Ennek az egyenletnek a -a21 szeresét adjuk hozzá a (2) egyenlethez, a -a31 szeresét a (3) egyenlethez, és a -a41 szeresét a (4) egyenlethez! Ekkor az új egyenlet rendszerben az első sort kivéve az x1 együtthatója mindenütt nulla lesz. Ezután végezzük el ugyan ezt a műveletet a második egyenlettel, de már az új együtthatókkal. Ekkor az egyenletrendszerben a második sor kivételével x2 együtthatói nullák. Folytatva az eljárást, azaz i=1,2,3,4-re az i-edik sort elosztjuk aii-vel, és minden j=/=i-re, a j-edik egyenlethez hozzáadva az i-edik egyenlet -aji szeresét, az i-edik sort kivéve xi együtthatói nullák lesznek. Mire elvégezzük a műveletet az összes sorral, minden egyes egyenletben pontosan egy darab együttható lesz nem nulla, hanem egy. Ekkor az egyenletrendszerünk így néz ki:
Ahol bi* az i-edik egyenlet jobb oldalán álló konstans tag a műveletsor végén. Látszik, hogy az egyenletrendszer megoldása: xi=bi* Ha nem csak egy konstans vektorunk van, akkor minden egyes műveletet el kell végezni minden egyes, az éppen adott sorban lévő konstans tagokkal. Inverz mátrix meghatározása Gauss--Jordan eliminációvalHa az egyenletrendszer megoldásakor használt együttható mátrix A, az ismeretlenek vektora x1 x2... xn és a konstans tagok vektora b1 b2... bn, akkor az egyenletrendszer az
alakba írható. ( Mind x, mind pedig b oszlop vektorok.) Ha nem egy ilyan egyenletünk van, hanem n darab, és A mindegyikben ugyanaz, akkor ezek felírhatók:
alakban. Tulajdonképpen nem teszünk mást, mint ugyan ahhoz az A együttható mátrixhoz, és különböző b konstans vektorokhoz keressük a megoldást. Ha az egymás mellé írt n darab x vektor helyett egy X mátrixot írunk, amelynek xij eleme éppen a j-edik x vektor i-edik komponense, és hasonló módon a b vektorok helyett egy B mátrixot írunk, akkor ismerve a mátrix szorzás szabályait könnyen belátható, hogy A X=B egyenlőség igaz, ha az egyes x vektorok megoldásai a megfelelő (5) egyenletnek. Ha a B mátrix speciális módon pont az egységmátrix, akkor az X mátrix az A mátrix inverze. Nem kell tennünk tehát semmi mást, mint követve a Gauss--Jordan eliminációt az n darab egyenletrendszert egyszerre megoldjuk, és az eljárás végén az egységmátrix helyén megjelenik az A mátrix inverze, A-1. A Gauss--Jordan elimináció során, akár egyenlet rendszer megoldására használjuk, akár inverz mátrix előállítására, előfordulhat, hogy a főátló beli elem, amellyel éppen osztani akarunk nulla. Ez nem jelenti feltétlenül azt, hogy az egyenletrendszer szinguláris. Például az
egyenletrendszer megoldása során, ha a második egyenlethez hozzáadjuk az első sor -1 szeresét, akkor az
egyenletrendszert kapjuk. (a11-gyel nem kellett osztanunk az első sort, mert a11=1 volt.) Itt a második sorban x2 együtthatója nulla, nem lehet vele osztani. Az egyenletrendszernek mégis van egyértelmű megoldása, azaz nem szinguláris. Mivel ez a probléma csak a főelem kiválasztás nélküli Gauss--Jordan eliminációnál fordulhat elő, és mivel a Gauss--Jordan elimináció főelem kiválasztás nélkül nagyon könnyen nam stabilissá válhat, ezért ennek a problémának a megoldásával nem érdemes foglalkozni. Főelem kiválasztásAz elimináció során főelemnek, (angolul pivotnak) nevezzük azt a főátlóbeli elemet amelyikkel éppen elosztjuk az aktuális sort. Az elimináció akkor lesz leginkább stabilis, ha ez az elem abszolút értékben a lehető legnagyobb. Ha a főelem kiválasztás során vagy csak az aktuális sorból, vagy csak az aktuális oszlopból választjuk ki a legnagyobb abszolút értékűt, akkor részleges főelem kiválasztást, ha pedig az egész mátrixban keressük a legnagyobb elemet, akkor teljes főelem kiválasztást végzünk. Ahhoz, hogy a kiválasztott elem a főátlóba kerüljön sor és oszlopcseréket hajtunk végre a mátrixban, arra azonban ügyelni kell, hogy ezeket a cseréket a program megfelelően adminisztrálja. Ellenkező esetben ugyanis oszlop cserék után előfordulhat, hogy a konstans tagok helyén megkapjuk ugyanaz ismeretlenek értékét, csak azt nem fogjuk tudni, hogy melyik érték melyik ismeretlenhez tartozik. Általában mind a részleges, mind pedíg a teljes főelem kiválasztásnál el lehet kerülni az oszlopok cseréjét. Mindkét fajta főelem kiválasztásnál csak olyan elemeket vehetünk figyelembe, amelyek olyan oszlopban és sorban állnak, amelyeket az elimináció még nem normált le, vagy nullázott ki. Furcsának tűnhet a részleges főelem kiválasztásnál, hogy a legnagyobb abszolút értékű elemet választjuk ki. Az egyenletrendszerünk megoldása ugyanis nem változik attól, ha az egyik egyenletünket egy konstans számmal megszorozzuk. Egy megfelelően nagy számot választva szorzónak bármelyik sor eleméből csinálhatunk főelemet. Emiatt néha alkalmazzák az ugynevezett implicit főelem kiválasztást, amelyik nem a legnagyobb elemet választja ki az adott oszlopból, hanem azt, amelyik a legnagyobb lenne akkor, ha minden sort végig osztottunk volna a legnagyobb elemével. A főelem kiválasztásnál megtehetjük, hogy nem cseréljük meg a sorok elemeit, hanem egy integer array/vel tartjuk nyilván a cseréket. Ez a módszer nyilván gyorsabb, mint a sorok fizikai megcserélése, de ha az eredeti mátrix helyén elő akarjuk állítani az inverz mátrixot, akkor nem tudjuk a sorok valódi megcserélését kikerülni. Gauss elimináció visszahelyettesítésselVáltoztassunk a Gauss--Jordan elimináción úgy, hogy az egyes sorokat csak a nagyob sorszámú sorokhoz adjuk hozzá, azaz ha éppen az i-edik sort adjuk hozzá a többihez, akkor minden j>i-re a j-edik sorhoz hozzáadjuk az i-edik sor -aji/aii szeresét. Ekkor az elimináció végére az A mátrix főátló alatti része lesz kinullázva.
Az eljárást eddíg a pontig hívják Gauss eliminációnak. Ezután jön a visszahelyettesítés. Az egyenletrendszer nem kapott olyan triviális formát, mint a Gauss Jordan elimináció után, de xn értékét azért könnyen megkapjuk: xn=bn/ann Ha pedig ismerjük xn-et, akkor xn-1:
Általánosságban azt mondhatjuk, hogy ha ismerjük xj értékét j=i+1,..,n-re, akkor
Az az eljárás, amelyik ennek a képletnek az alapján határozza meg az ismeretleneket visszahelyettesítés. Ennek a módszernek az az előnye a Gauss--Jordan eliminációval szemben, hogy kevesebb művelet elvégzésére van szükség. Míg a Gauss--Jordan elimináció legbelső ciklusa, amely egy szorzást és egy kivonást tartalmaz n3-szor hajtódik végre, addig itt, mivel nem tisztítjuk ki az egész mátrixot, hanem csak a felét, 1/3 n3 -szor. Az a hátrány, hogy a konstans vektorokat mind ismernünk kell az elimináció megkezdésekor itt is jelen van. A Gauss--Jordan elimináció programjaA programot lépésenként fejlesztjük ki. Nézzük meg az első változatot! 001 procedure gaussj1( var a:glnpbynp; n:integer; 002 var b:glnpbymp; m:integer ); 003 (* 004 TYPE glnpbynp=ARRAY[1..np,1..np] of REAL; 005 glnpbymp=ARRAY[1..np,1..mp] of REAL; 006 *) 007 VAR i,ll,l:integer; 008 inv,dum:real; 009 BEGIN 010 for i:=1 to n do 011 BEGIN 012 inv := 1.0 / a[i,i]; 013 for l:=1 to n do 014 a[i,l] := a[i,l] * inv; 015 for l:=1 to m do 016 b[i,l] := b[i,l] * inv; 017 018 for ll:=1 to n do 019 if lli then 020 BEGIN 021 dum := a[ll,i]; 022 for l:=1 to n do 023 a[ll,l] := a[ll,l] - a[i,l] * dum; 024 for l:=1 to m do 025 b[ll,l]:= b[ll,l] - b[i,l] * dum ; 026 END 027 END 028 END; Ez a program nem végez főelem kiválasztást, és nem is állítja elő az a mátrix inverzét, csak kinullázza a főátló alatti és fölötti elemeket. A legkülső ciklus I végigmegy minden soron. A cikluson belül először elosztja az éppen aktuális sor minden elemét (a konstans tagokat is) a főátló beli elemmel. Ezzel tulajdonképpen ezt a sort normáltuk, és a főátlóba 1.0 került. A következő ciklus ll végigmegy minden soron, kivéve az i-edik sort, és minden sorhoz hozzáadja az i-edik sor annyiszorosát, hogy az i -edik oszlopba nulla kerüljön. Ez a ciklus kinullázza az éppen aktuális sor főátlójában lévő elem alatti és feletti elemeket. Az i -edik sorból természetesen nem szabad kivonni önnmagát, hiszen ezzel információt veszítenénk, kinulláznánk az egész sort, a konstans tagokat is, és nem tudnánk megoldani az egyenlet rendszert. Az eljárás lefutása után az a mátrix helyén megkapjuk az egységmátrixot, a b elemeiben pedig az ismeretlenek értékeit. Hogyan állítsuk elő az A mátrix inverzét? Ha a B mátrix mérete n*n, és eredeti tartalma b= 001 procedure gaussj2( var a:glnpbynp; n:integer ); 002 (* 003 TYPE glnpbynp=ARRAY[1..np,1..np] of REAL; 004 *) 005 VAR i,ll,l:integer; 006 inv,dum:real; 007 BEGIN 008 for i:=1 to n do 009 BEGIN 010 inv := 1.0 / a[i,i]; 011 a[i,i] := 1.0; 012 for l:=1 to n do 013 a[i,l] := a[i,l] * inv; 014 015 for ll:=1 to n do 016 if lli then 017 BEGIN 018 dum := a[ll,i]; 019 a[ll,i] := 0.0; 020 for l:=1 to n do 021 a[ll,l] := a[ll,l] - a[i,l] * dum; 022 END 023 END 024 END; Az eljárás új változatában nem kell tehát törődnünk a b mátrixszal külön, hiszen a megfelelő részeit mind az a, mind pedig a b mátrixnak az a array-ben tároljuk. A legkülső ciklus i-edik lefutásakor, miután a főátlóbeli elemet kivettük a mátrixból, és a reciprokát az inv változóban helyeztük, a főátlóba egyet írunk, így a sor végig osztása után ide pont a főátló beli elem reciproka kerül, vagyis az, ami az inverz mátrix kialakítása során a bi pillanatnyi értéke kell, hogy legyen. Az ll ciklusban, amikor a[ll,i]-t kinulláztuk, éppen az inverz mátrix [ll,i] pozícióban levő értékét állította elő az a array [ll,i] elemében a program. Amikor l\le i, akkor mindkét l ciklus az inverz mátrixot állítja elő, amikor pedig i<l, akkor az eliminációt végzi. Természetesen miközben előállítjuk az a array-ben az inverz mátrixot a b array-ben az egyenletrendszer megoldásait is létrehozhatjuk. 001 procedure gaussj3( var a:glnpbynp; 002 n:integer; 003 var b:glnpbymp; 004 m:integer ); 005 (* 006 TYPE glnpbynp=ARRAY[1..np,1..np] of REAL; 007 glnpbymp=ARRAY[1..np,1..mp] of REAL; 008 *) 009 VAR i,ll,l:integer; 010 inv,dum:real; 011 BEGIN 012 for i:=1 to n do 013 BEGIN 014 inv := 1.0 / a[i,i]; 015 a[i,i] := 1.0; 016 for l:=1 to n do 017 a[i,l] := a[i,l] * inv; 018 for l:=1 to m do 019 b[i,l] := b[i,l] * inv; 020 021 for ll:=1 to n do 022 if lli then 023 BEGIN 024 dum := a[ll,i]; 025 a[ll,i] := 0.0; 026 for l:=1 to n do 027 a[ll,l] := a[ll,l] - a[i,l] * dum; 028 for l:=1 to m do 029 b[ll,l]:= b[ll,l] - b[i,l] * dum ; 030 END 031 END 032 END; Meg is lennénk az eliminációval, ha ezek a programok nem lennének veszélyesek. A veszély abban rejlik, hogy nem végeztünk főelem kiválasztást, és ezért a halmozódó hibák miatt nagyon könnyen a valódi megoldástól távoli eső számokat kaphatunk. Mielőtt azonban nekilátnánk a főelem kiválasztásnak próbáljuk meg végig gondolni, hogy vajon működnének-e a programjaink, ha a legkülső ciklusban i értéke nem egytől menne egyesével n-ig, hanem az 1,2, ... ,n számokat valamely más permutáció szerint venné sorra. Természetesen a programok ekkor is ugyanazt az eredményt adják. Meg kell még jegyezni azt is, hogy egyik program sem vizsgálja, hogy a főátlóbeli elem, amivel a sort végigosztjuk nem nulla-e. Ha ez az elem nulla, az még nem jelenti azt, hogy a mátrix szinguláris. (Ahogy azt korábban már egy példán keresztül bemutattuk.) A főelem kiválasztást alkalmazó programokat is lépésenként fejlesztjük ki. Készítsük el először azt a programot, amelyik részleges főelem kiválasztást végez, és nem állítja elő az inverz mátrixot! A főelemet az éppen aktuális elem oszlopában keressük. Azért az oszlopban, és nem a sorban, mert mert így ahhoz, hogy a kiválasztott elem arra a főátló beli helyre kerüljön ahol éppen tartunk, sorokat és nem oszlopokat kell megcserélni. Ha oszlopokat cserélnénk meg, akkor nyilván kellene tartani a cseréket, egy egyenletrendszer megoldását azonban az egyenletek sorrendje nem befolyásolja. A sorcseréket természetesen a B mátrixban is végre kell hajtani, hiszen a konstans tagok is hozzá tartoznak az egyenletekhez. Ha az oszlopban a legnagyobb abszolút értékű elem nulla, akkor minden elem nulla, és ez azt jelenti, hogy az együttható mátrix szinguláris, az egyenletrendszernek nincs egyértelmű megoldása. Az egyes sorok cseréje után a program hátralevő része ugyanaz, mint a főelem kiválasztás nélküli esetnél. 001 procedure gaussj4( var a:glnpbynp; n:integer; 002 var b:glnpbymp; m:integer ); 003 (* 004 TYPE glnpbynp=ARRAY[1..np,1..np] of REAL; 005 glnpbymp=ARRAY[1..np,1..mp] of REAL; 006 *) 007 VAR i,ll,l,row:integer; 008 inv,dum,big:real; 009 BEGIN 010 for i:=1 to n do 011 BEGIN 012 row := 0; 013 big := 0.0; 014 for l:=1 to n do 015 if big <= abs(a[l,i]) then 016 BEGIN 017 row = l; 018 big= abs(a[l,i]); 019 END; 020 if big="0.0" then 021 BEGIN 022 writeln('The matrix is singular.'); 023 halt; 024 END; 025 if row> i then 026 BEGIN 027 for l:=1 to n do 028 BEGIN 029 dum := a[i,l]; 030 a[i,l] := a[row,l]; 031 a[row,l] := dum; 032 END; 033 for l:=1 to m do 034 BEGIN 035 dum := b[i,l]; 036 b[i,l] := b[row,l]; 037 b[row,l] := dum; 038 END; 039 END; 040 inv := 1.0 / a[i,i]; 041 for l:=1 to n do 042 a[i,l] := a[i,l] * inv; 043 for l:=1 to m do 044 b[i,l] := b[i,l] * inv; 045 for ll:=1 to n do 046 if lli then 047 BEGIN 048 dum := a[ll,i]; 049 for l:=1 to n do 050 a[ll,l] := a[ll,l] - a[i,l] * dum; 051 for l:=1 to m do 052 b[ll,l]:= b[ll,l] - b[i,l] * dum ; 053 END; 054 END; 055 END; Mi a helyzet akkor, ha az inverz mátrixot is elő akarjuk állítani? Amikor az A X=I egyenletet oldjuk meg, akkor minden egyes sorművelet amit az A és az I mátrixon egyidőben végzünk megfelel valamilyen P mátrixszal való szorzásnak. Az egész elimináció egyenlete:
ahol Pi az i-edik sorműveletnek megfelelő mátrix. A műveletek során a mátrix szorzásokat jobbról balra haladva végezzük el, de mivel a mátrix szorzás asszociatív semmi sem tiltja, hogy gondolatban most balról jobbra végezzük el, és a
jelölést használva:
A főelem kiválasztás sorcseréi ugyan olyan mátrix műveletek, mint az elimináció többi sorművelete, így igaznak kell lennie, hogy minden további sorművelet nélkül
Van azonban egy kis probléma. Nevezetesen az, hogy mi nem állítjuk elő az I mátrixot a kezdet kezdetén, hanem csak menet közben az elimináció során. Ha pedig most bővítve a gaussj4 programunkat a 040 sor utánra beszúrjuk az a[i,i]:=1.0; sort, s a 048 sor utánra pedig az a[ll,i]:=0.0; sort, akkor az egységmátrix i-edik oszlopának az elóállítása már a sorok cseréje után történt meg, és így az ebben az oszlopban álló elemet nem cseréltük meg, és mivel minden j>i j-edik oszlopot is később állítunk elő, a két sor minden j>=i j-edik elemét nem cseréltük meg. Szerencsére ezeknem az elemeknek a zöme nulla, csupán az i-edik sor i-edik eleme, és ha az a sor amelyikkel az i-edikel megcseréltük a j-edik, akkor ennek a j-edik eleme nem nulla, hanem 1. Így, hogyha az i-edik, és j-edik sort megcseréljük, majd pedig az i-edik és j-edik oszlopot, akkor pont azt a mátrixot kapjuk, amelyikben az i-edik és j-edik sor első i-1 eleme lett megcserélve. Ezzel tehát, hogy az 1-est a főátlóba, és a nullákat az éppen aktuális oszlop helyére a sor csere után írtuk be tulajdonképpen egy oszlop cserét hajtottunk végre az egységmátrixból az inverz mátrixszá átlakuló mátrix azon részén amelyik még érintetlenül tartalmazza az egységmátrix elemeit. Mivel azonban ezt az oszlopcserét nem hajtottuk végre az egyenlet bal oldalán, ezért még egyszer végre kell hajtanunk, hogy a hatását megszüntessük. Az egyes oszlopcserék a mátrix egyenletbe az egységmátrix jobbról való Q mátrixszokkal való szorzásával írhatók be.
|