Numerikus AnalízisLU 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:
Ezt a felbontást felhasználhatjuk egy egyenletrendszer megoldására:
és i=2,3, ... ,n-re
A második pedig visszafelé haladó behelyettesítéssel:
és i=n-1,n-2, ... ,2,1-re
Ez a megoldás n2-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 n2 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 n2+n darab van. (Az
(Figyeljük meg, hogy a (6) és (7) képletek csak a 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 Az LU dekompozíció programjaAz 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 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: ![]() |