LU-Zerlegung

Die LU-Zerlegung, auch LR- oder Dreieckszerlegung ist ein Algorithmus zur Programmierung des Gauss-Jordan-Verfahrens. Hier wier die Idee anhand einer  Untersuchung im Falle einer 3 x 3-Matrix gezeigt. Die Methode läst sich dann mühelos auf den Fall einer beliebigen n x n-Matrix erweitern.

Bedeutung von LU: L  Lower, U Upper.                                                  © Rolf Wirz

Hier verwendete Programmversion: Mathematica 5.2

1. Zur Theorie der LU-Zerlegung

Das Beispiel einer abstrakten 3 x 3-Matrix:

Konstruktion von U  (dabei werden die Elementarsubstitutionen werden in Matrizenoperationen abgebildet):

Sei a1 nicht 0 (sowie weiter unten auch die Nenner nie 0 - ansonst eine Zeilenvertauschung vorgenommen werden müsste...):

Remove["Global`*"]

A30 = {{a1,b1,c1},{a2,b2,c2},{a3,b3,c3}}; 
K31 = {{0,0,0},{0,0,0},{A30[[3]][[1]]/A30[[1]][[1]],0,0}};
Map[MatrixForm,{A30, K31}]

H13 = IdentityMatrix[3]- K31; H13 // MatrixForm

U1 = H13.A30//Simplify; U1 //MatrixForm

K21 = {{0,0,0},{U1[[2]][[1]]/U1[[1]][[1]],0,0},{0,0,0}}; K32 // MatrixForm

H12 = IdentityMatrix[3]-K21; H12 // MatrixForm

U2 = H12.U1//Simplify; U2 //MatrixForm

Sei -((a2 b1)/a1) + b2 nicht 0:

K32 = {{0,0,0},{0,0,0},{0,U2[[3]][[2]]/U2[[2]][[2]],0}}; K32 // MatrixForm

H23 = IdentityMatrix[3]-K32; H23 // MatrixForm

U3 = H23.U2//Simplify; U3 //MatrixForm

In einem Schritt:

U3 = H23.H12.H13.A30//Simplify; U3 //MatrixForm

Uresult = U3;

Konstruktion von L mittels inverser Matrix:

L U = A ==> L = A InverseU. Daher ist die Inverse der Dreiecksmatrix U zu bestimmen. 
Ansatz:  U * InverseU = E. Die Berechnung von InverseU ist hier einfach. Die Sukzessive Berechnung der Elemente von InverseU zeigt rasch, dass InverseU ebenfalls eine obere Dreiecksmatrix sein muss.

U = {{u11,u12,u13},{0,u22,u23},{0,0,u33}}; U // MatrixForm

InvU = {{v11,v12,v13},{0,v22,v23},{0,0,v33}}; InvU // MatrixForm

U.InvU// MatrixForm

InvU.U// MatrixForm

solv = Solve[U.InvU == IdentityMatrix[3],{v11,v12,v13,v22,v23,v33}] // Flatten

Hier hat man ein lineares Gleichungssystem mit 6 Unbekannten, das man nach dem Rückwärtseinsetzungsverfahren rasch und problemlos lösen kann.

InvU = InvU /.solv; InvU // MatrixForm

Kontrolle:

Inverse[U] // MatrixForm

Ersetzung der künstlichen Koeffizienten in U durch die von A30:

UFlat = Flatten[U];
U3Flat = Flatten[U3];
rul = Table[UFlat[[k]]->U3Flat[[k]],{k,1,Length[UFlat]}]

L = A30.InvU; Lresult = L/.rul //Simplify; Lresult // MatrixForm

Kontrolle:

Lresult.Uresult //Simplify // MatrixForm

Simplify[Lresult.Uresult] == A30

Anwendung für die Lösung eines Gleichungssystems:

Sei    A.x = (L.U).x = L.(U.x) = L.y = b .
Löse zuerst   L.y = b   ==>  y wird bekannt (einfache Lösungsmöglichkeiten, da L Dreiecksmatrix ist).
Löse anschliessend   U.x = y   ==>  x wird bekannt (einfache Lösungsmöglichkeiten, da U ebenfalls Dreiecksmatrix ist).

0. Materialbereitstellung

Remove["Global`*"]

A3 = {{1,2,3},{1,-1,0},{-2,1,5}}; 
B3 = {{3,5,1},{4,7,9},{3,2,6}}; 
C3 = {{4,1,3},{5,6,5},{5,8,8}};
Map[MatrixForm,{A3,B3,C3}]

X31 = {{x11},{x21},{x31}}; 
Y31 = {{y11},{y21},{y31}}; 
B31 = {{1},{2},{3}}; 
Map[MatrixForm,{X31,Y31,b31}]

1. LU-Zerlegung für 3 x 3 Matrizen als Modul

(Einfach lesbarer Modul ohne kompakte Programmierung, wiederholt anwendbar!)

modulLU[{{a1_,b1_,c1_},{a2_,b2_,c2_},{a3_,b3_,c3_}}]:=
Modul[{},
A30 = {{a1,b1,c1},{a2,b2,c2},{a3,b3,c3}}; 
K31 = {{0,0,0},{0,0,0},{A30[[3]][[1]]/A30[[1]][[1]],0,0}};
Map[MatrixForm,{A30, K31}];
H13 = IdentityMatrix[3]- K31;
U1 = H13.A30//Simplify;
K21 = {{0,0,0},{U1[[2]][[1]]/U1[[1]][[1]],0,0},{0,0,0}}; 
H12 = IdentityMatrix[3]-K21; 
U2 = H12.U1//Simplify;
K32 = {{0,0,0},{0,0,0},{0,U2[[3]][[2]]/U2[[2]][[2]],0}};
H23 = IdentityMatrix[3]-K32;
U3 = H23.U2//Simplify;
Uresult = U3;
U = {{u11,u12,u13},{0,u22,u23},{0,0,u33}};
InvU = {{v11,v12,v13},{0,v22,v23},{0,0,v33}}; 
solv = Solve[U.InvU == IdentityMatrix[3],{v11,v12,v13,v22,v23,v33}] // Flatten;
InvU = InvU /.solv; 
UFlat = Flatten[U];
U3Flat = Flatten[U3];
rul = Table[UFlat[[k]]->U3Flat[[k]],{k,1,Length[UFlat]}];
L = A30.InvU; Lresult = L/.rul //Simplify;
Print["Eingabematrix = ",A30//MatrixForm];
solv1:= Solve[Evaluate[Flatten[Lresult.Y31]  == Flatten[B31] ],{y11,y21,y31}]//Flatten;
YN31=Y31/.solv1;
solv2:= Solve[Evaluate[Flatten[Uresult.X31]  == Flatten[YN31] ],{x11,x21,x31}]//Flatten;
XN31=X31/.solv2;
Print["U = ",Uresult//MatrixForm];
Print["L = ",Lresult//MatrixForm];
Print["Kontrolle: L U = ",Lresult.Uresult//MatrixForm];
Print["Y21 = ",YN31//MatrixForm];
Print["X21 = ",XN31//MatrixForm];]
modulLU[A3];

Kontrolle:

Solve[A3.Flatten[X31]==Flatten[B31],Flatten[X31]]

Neue Anwendung des Moduls:

modulLU[B3];

Kontrolle:

Solve[B3.Flatten[X31]==Flatten[B31],Flatten[X31]]

Neue Anwendung des Moduls:

modulLU[C3];

Solve[C3.Flatten[X31]==Flatten[B31],Flatten[X31]]