Cholesky Zerlegung in Matlab

Question

Aufgabe:

Ich soll in Matlab die Cholesky-Zerlegung programmieren. Bisher bin ich so weit:

function [L]= cholesky(A)

if ~issymmetric(A)
   % Das Zeichen ~ ist die Negation eines logischen Ausdruckes in Matlab 
   disp('Die Matrix ist nicht symmetrisch.')
   % Hier muss ein Warnhinweis gegeben werden!
   return
   % return beendet hier die Funktion
end

[n]=size(A,1); 
%size gibt hier die Anzahl der Zeilen an (Zeilen wegen der 1, bei
%Spalten 2 aber ist hier egal)

L=zeros(n);

for k=1:n
   L(k,k)=(A(k,k)-sum(L(k,1:k-1).^2))^(1/2);
   for i=k+1:n
       L(i,k)=(A(i,k)-sum(L(i,1:k-1).*L(k-1,1:k-1)))/L(k,k)
   end
end

if ~isreal(L)
   disp('Achtung: Die Matrix ist nicht positiv definit!')
   % Wenn die Matrix komplexe Werte enthält, war die Cholesky Zerlegung
   % eigentlich nicht durchführbar und somit ist die Matrix nicht positiv
   % definit
end


end

Problem/Ansatz:

In der markierten Zeile wird mir beim Gleichheitszeichen ein Fehler angezeigt. Weiß jemand weshalb? Es ist meine erste Programmierhausaufgabe, ich habe also noch gar keine Erfahrung mit Matlab.

Danke schon mal =)

Answer 1

Wenn ich die Matrix (sagen wir R(4,4))

\(\scriptsize L \, :=  \, \left(\begin{array}{rrrr}l11&0&0&0\\l21&l22&0&0\\l31&l32&l33&0\\l41&l42&l43&l44\\\end{array}\right)\)

mit der Transponierten multiplziere soll A {a_ij} rauskommen, also L L^T - A =0

\(\scriptsize LLT=\left(\begin{array}{rrrr}l11^{2} - a11&l11 \; l21 - a12&l11 \; l31 - a13&l11 \; l41 - a14\\l11 \; l21 - a21&l21^{2} + l22^{2} - a22&l21 \; l31 + l22 \; l32 - a23&l21 \; l41 + l22 \; l42 - a24\\l11 \; l31 - a31&l21 \; l31 + l22 \; l32 - a32&l31^{2} + l32^{2} + l33^{2} - a33a34&l31 \; l41 + l32 \; l42 + l33 \; l43\\l11 \; l41 - a41&l21 \; l41 + l22 \; l42 - a42&l31 \; l41 + l32 \; l42 + l33 \; l43 - a43&l41^{2} + l42^{2} + l43^{2} + l44^{2} - a44\\\end{array}\right)=0\)

da würde ich spaltenweise runter rechnen...

Vergleich das mal mit Deinem Code - könnte das sein das da ein Indexdreher drin ist?

Nachtrag: da entsteht irgendwann mal ...L(1,1:0)*L(0,1:0)...wie wird das interpretiert?

Comments

Zum Nachtrag: Das habe ich jetzt auch überprüft, aber zwei Zeilen davor funktioniert das auch. Es wird einfach eine 0 addiert.

Wenn ich die innere for-Schleife raus lasse, dann läuft das Programm durch. Also es liegt auch sicher an der Zeile. Der Fehlercode, den ich erhalte, ist folgender:

Index in position 1 is invalid. Array indices must be positive integers or logical values.

Error in cholesky (line 20)
      L(i,k)=(A(i,k)-(sum(L(i,1:k-1).*L(k-1,1:k-1))))/L(k,k);

Ich weiß aber nicht, was das bedeutet.

Die Indizes müssten alle stimmen..

---

>Index in position 1 is invalid. Array indices must be positive integers or logical values.

Es liegt aber an der Indexierung und da ist L(k-1,1:k-1) ein guter Kandidat, k=1 macht L(0,...)

>Index in position 1..

Ich hab kein MatLab nur was cloneartiges da schreib ich

Bei den L_jj,L_ij könnte die Summe bis k=1..j-1 gehen wenn es keine Indexprobleme macht...

Man könnte auch L₁₁ vorbelegen und L_{J+1 J+1} am Ende der Schleife definieren?

---

Genau an der Stelle hat es gehakt, danke für die Hilfe!

Answer 2

Hallo,

wenn du L(k-1,1:k-1) in L(k,1:k-1) umänderst sollte der Code funktionieren.

Dann lautet die ganze Zeile:

for i=(k+1):n
      L(i,k)= (A(i,k)-sum(L(i,1:k-1).*L(k,1:k-1)))/L(k,k);

Comments

Das war es, danke!