Numerikus Analízis

LU dekompozíció

Tegyük fel, hogy az A mátrixot sikerült L U = A alakban felírni, ahol L alsó háromszög mátrix, azaz csak a főátlóban és az alatt tartalmaz nem nulla elemeket, U pedig felső háromszög mátrix, azaz csak a főátlóban és a felett tartalmaz nem nulla elemeket. Egy A 4*4-es mátrixra ez így néz ki:

(1)

Ezt a felbontást felhasználhatjuk egy egyenletrendszer megoldására: úgy, hogy először megoldjuk az egyenlet rendszert, majd pedig az egyenletrendszert. Mi az értelme annak, hogy egy egyenletrendszer helyett kettőt oldunk meg? Az, hogy mivel mind az L, mind pedig az U háromszög mátrixok, ezért ennek akét egyenletrendszernek a megoldása triviális. Az első egyenlet megoldható előre haladó behelyettesítéssel:

(2)

és i=2,3, ... ,n-re

(3)

A második pedig visszafelé haladó behelyettesítéssel:

(4)

és i=n-1,n-2, ... ,2,1-re

(5)

Ez a megoldás n²-tel arányos mennyiségű számítást igényel, és valahányszor egy újabb b vektorral kívánjuk megoldani az egyenletet csak a (2) ... (5) képleteket kell alkalmazni. Most már csak az a kérdés, hogy egy adott A mátrixhoz hogyan keressük meg a megfelelő L és U mátrixokat. Ha az (1) mátrix egyenletben elvégezzük a beszorzást, akkor n² darab skalár egyenletet kapunk. Írjunk fel ebből egyet!

A három pont itt nem azt jelenti, hogy az összegzést vég nélkül kell folytatni, hanem azt, hogy i és j értékétől függően másképpen kell befejezni. Ha ugyanis i<j, akkor:

ugyanis minden \beta_{kj}=0 ha k>j. Ha i=j:

és ha i>j, akkor:

Ez pontosan n darab egyenletet jelent, ismeretlenünk pedig n²+n darab van. (Az és értékek teljesen kitöltik az n*n-es mátrixot, de a főátló kétszer szerepel, innen van a +n.) Hogy meglegyen a kellő számú egyenletünk, írjunk fel még n darabot! Legyen i=1,2, ... ,n-re Most tehát van n²+n darab egyenletünk, és ugyanennyi ismeretlenünk. Persze elég furcsának tűnhet, hogy egy n ismeretlenes egyenletrendszer helyett egy n²+n ismeretlenes, tehát jóval nagyobb egyenletrendszert akarunk megoldani. A magyarázat az, hogy bár ez az egyenletrendszer jóval nagyobb, mégis könnyen megoldható, mert speciális. A megoldás Crout algoritmusa:

• Legyen és . (Mivel i=j=1-re )
• Ezek után minden j=2,3, ... ,n értékre:
• Először i=1,2, ... ,j-re keressük meg a -t az egyenleteink alapján!

(6)

• (Ha i=1, akkor definiáljuk a értékét nullának.)
• Másodszor pedig i=j+1,j+2, ... ,n-re keressük meg _ij-t:

 

(7)

(Figyeljük meg, hogy a (6) és (7) képletek csak a való osztásban térnek el egymástól.) Könnyen látható, hogy követve az eljárást a következő és értékek kiszámításához a korábbi és értékek rendelkezésre állnak. Az is látszik, hogy amikor egy _ij-t vagy -t kiszámítottunk, akkor ahhoz felhasználjuk a_ij értékét, de arra a továbbiakban nincsen szükség a dekompozícióhoz, ezért az eredeti mátrixban tárolhatjuk a kiszámított és értékeket. (Ez alól kivételek az _ii értékek, ezeket azonban egyáltalán nem tároljuk, tudván, hogy ezek értéke 1.)

A főelem kiválasztást a stabilitás érdekében itt is alkalmaznunk kell, hatékonyan azonban csak részleges főelem kiválasztást tudunk megvalósítani. (A részleges főelem kiválasztás a stabilitáshoz legtöbb esetben ugyan olyan jó, mint a teljes főelem kiválasztás.) A részleges főelem kiválasztás alkalmazása itt azt jelenti, hogy nem magát az A mátrixot bontjuk fel L és U mátrixok szorzatára, hanem egy, az A mátrixból sorcserékkel előállított mátrixét. Természetesen ezeket a sorcseréket a felbontás későbbi felhasználása érdekében nyilván kell tartanunk. Itt sajnos nem tehetjük meg, hogy az eljárás végén helyrehozzuk a cseréket. (Ugyanis, míg egy A mátrixnak mindig van A^-1 inverze, ha nem szinguláris, LU felbontása nem biztos, hogy van.) A főelem kiválasztás megoldja azt a problémát is, ami akkor áll elő, ha az eredeti A mátrixnak nincs LU felbontása, bár nem szinguláris. A részleges főelem kiválasztást úgy végezzük el, hogy az a elem, amivel osztunk az _ij-t előállító képletben, a lehető legnagyobb legyen. Ehhez először kiszámítjuk a értékeket, majd az _ij értékek -szeresét. (Nem osztunk -vel az (7) képletben.) Ezután végezzük el a főelem kiválasztást. Kiválasztjuk a legnagyobb elemet a j-edik oszlopban az i>=j elemek közül, és az ilyen módon kiválasztott i-edik sort megcseréljük a j-edikkel. Ez a sorcsere a mátrix l=1,2, ... ,j-1 oszlopában a már felbontott rész _jl elemeit cseréli meg az _il elemekkel, és ez ugyanaz, mintha az eredeti mátrix elemeit cseréltük volna meg még a felbontás megkezdése előtt. A mátrix l=j+1, ... ,n oszlopaiban az eredeti mátrix elemei még sértetlenek. A j-edik oszlop megcserélésével szintén ugyanazt a hatást érjük el, mintha még a felbontás megkezdése előtt cseréltük volna meg a két sort. A két sor megcserélése után pedig a j-edik oszlop j+1, ... ,n elemeit elosztjuk a kiválasztott főelemmel.

Az LU dekompozíció programja

Az itt közölt program implicit főelem kiválasztást végez. Ezért még mielőtt hozzálátna a dekompozícióhoz, minden sorból kikeresi a legnagyobb abszolút értékű elemet, és az abszolút érték reciprokát elhelyezi a vv vektor megfelelő elemében. Ezután lát neki a dekompoz?ciónak. Végig megy minden oszlopon. Minden oszlopban először kiszámítja a értéket i=1,2, ... ,j-1-re, majd az új _ij értékeket i=j, ... ,n-re , de nem végzi el a -vel való osztást, hanem előbb kiválasztja az oszlop alsó részéből a legnagyobb abszolút értékűt (szorozva mindig a sor normájával), és megcseréli a sorokat. imax tartalmazza annak a sornak a számát, amelyet meg kell cserélni a j-edikkel. Ezután a vv vektor elemeit is meg kell cserélni, mivel azonban a j-edikre többet nem lesz szükség, ezért a csere műveletének a kétharmada elhagyható. Ezután a program elvégzi a kiválasztott főelemmel való osztást. Az ludcmp eljárásnak be, és kimenő paramétere az a mátrix. Ennek a méretét határozza meg az n egész változó. Az indx egészeket tartalmazó vektor kimenő változó, míg d értéke az eljárás lefutása után +1, vagy -1 attól függően, hogy páros, vagy páratlan számú sorcserét végeztünk.

PROCEDURE ludcmp(VAR a: glnpbynp; n: integer; VAR indx: glindx; VAR d: real);
CONST tiny=1.0e-20;
VAR k,j,imax,i: integer;
sum,dum,big: real; vv: glnarray;
BEGIN d := 1.0; FOR i := 1 to n DO
 BEGIN
 big := 0.0;
 FOR j := 1 to n DO
 IF abs(a[i,j]) > big THEN big := abs(a[i,j]);
 IF big = 0.0 THEN
 BEGIN
 writeln('pause in LUDCMP - singular matrix');
 halt;
 END; vv[i] := 1.0/big;
 END;
 FOR j := 1 to n DO
 BEGIN
 IF j > 1 THEN
 BEGIN
 FOR i := 1 to j-1 DO
 BEGIN
 sum := a[i,j];
 IF i > 1 THEN
 BEGIN
 FOR k := 1 to i-1 DO
 BEGIN
 sum := sum-a[i,k]*a[k,j];
 END;
 a[i,j] := sum;
 END;
 END;
 END;
 big := 0.0;
 FOR i := j to n DO
 BEGIN
 sum := a[i,j];
 IF j > 1 THEN
 BEGIN
 FOR k := 1 to j-1 DO
 BEGIN
 sum := sum-a[i,k]*a[k,j];
 END;
 a[i,j] := sum;
 END;
 dum := vv[i]*abs(sum);
 IF dum > big THEN
 BEGIN
 big := dum;
 imax := i;
 END;
 END;
 IF j imax THEN
 BEGIN
 FOR k := 1 to n DO
 BEGIN
 dum := a[imax,k];
 a[imax,k] := a[j,k]
 a[j,k] := dum;
 END;
 d := -d;
 vv[imax] := vv[j];
 END;
 indx[j] := imax;
 IF j n THEN
 BEGIN
 IF a[j,j] = 0.0 THEN a[j,j] := tiny;
 dum := 1.0/a[j,j];
 FOR i := j+1 to n DO
 BEGIN
 a[i,j] := a[i,j]*dum;
 END;
 END;
 END;
 IF a[n,n] = 0.0 THEN a[n,n] := tiny
END; 

Az indx vektort kell, hogy használja az a program, amely a felbontott mátrix segítségével megoldja az egyenletrendszert. A d változó az a determinánsának a kiszámításához lehet szükséges. Ha az A mátrixot felbontottuk az L és U mátrixok szorzatára, akkor az eredeti egyenletrendszert nagyon egyszerűen, a (2) ... (5) képletek segítségével lehet megoldani. Ezt végzi el az lubksb eljárás.

PROCEDURE lubksb( a: glnpbynp; n: integer; indx: glindx; VAR b: glnarray);
VAR j,ip,ii,i: integer; sum: real;
BEGIN
 ii := 0;
 FOR i := 1 to n DO
 BEGIN
 ip := indx[i];
 sum := b[ip];
 b[ip] := b[i];
 IF ii 0 THEN
 FOR j := ii to i-1 DO
 sum := sum-a[i,j]*b[j]
 ELSE IF sum 0.0 THEN ii := i; b[i] := sum;
 END;
 FOR i := n DOWNTO 1 DO
 BEGIN
 sum := b[i];
 IF i<n THEN
 FOR j i+1 to n DO
 sum := sum-a[i,j]*b[j];
 b[i] := sum/a[i,i];
 END;
END;

Ha ez a két eljárás a rendelkezésünkre áll, akkor egy egyenletrendszer megoldása két sorból áll:

ludcmp(A,N,indx,glindx,d);
lubksb(A,N,indx,b)

Ha pedig a későbbiek folyamán ugyanezekkel az együtthatókkal, de egy c konstans vektorral kívánjuk megoldani az egyenletet, akkor már csak a második eljárást kell meghívnunk:

lubksb(A,N,indx,c)

(Az elsőt nem is szabad!) Hogyan számítsuk ki a mátrix determinánsát, ha ismert az LU felbontás? A determinánsra igaz, hogy tetszőleges, de azonos méretű négyzetes mátrixok esetén: