2015-03-30 131 views
1

我試圖與主要基於LU decomposition with partial pivoting MatlabMATLAB LU分解的部分迴轉

function [L,U,P] = lup(A) 
n = length(A); 
L = eye(n); 
U = zeros(n); 
P = eye(n); 

for k=1:n-1 
% find the entry in the left column with the largest abs value (pivot) 
[~,r] = max(abs(A(k:end,k))); 
r = n-(n-k+1)+r; 

A([k r],:) = A([r k],:); 
P([k r],:) = P([r k],:); 
L([k r],:) = L([r k],:); 

% from the pivot down divide by the pivot 
L(k+1:n,k) = A(k+1:n,k)/A(k,k); 

U(k,1:n) = A(k,1:n); 
A(k+1:n,1:n) = A(k+1:n,1:n) - L(k+1:n,k)*A(k,1:n); 

end 
U(:,end) = A(:,end); 

end 

這似乎對於大多數矩陣(等於MATLAB魯功能)工作,我的路分解工作,但下面的矩陣似乎以產生不同的結果:

A = [ 
3 -7 -2  2 
-3  5  1  0 
6 -4  0 -5 
-9  5 -5 12 
]; 

我只是不知道什麼是錯的。它似乎在鏈接提到的矩陣上工作正常

回答

2

你是非常接近。我換了三線總

for k=1:n-1成爲for k=1:n我們不這樣做的-1,因爲我們也希望得到L(n,n)=u(n,n)/u(n,n)=1用你的方法我們離開了這一點

L(k+1:n,k) = A(k+1:n,k)/A(k,k);成爲L(k:n,k) = A(k:n,k)/A(k,k);因爲你離開了L(k,k)=A(k,k)/A(k,k)=1

因爲k+1改變我們不需要先從對於L單位矩陣,因爲我們現在正在再現1點的對角線上,以便L=eyes(n);成爲L=zeros(n);

和完成代碼

function [L,U,P] = lup(A) 
% lup factorization with partial pivoting 
% [L,U,P] = lup(A) returns unit lower triangular matrix L, upper 
% triangular matrix U, and permutation matrix P so that P*A = L*U. 
    n = length(A); 
    L = zeros(n); 
    U = zeros(n); 
    P = eye(n); 


    for k=1:n 
     % find the entry in the left column with the largest abs value (pivot) 
     [~,r] = max(abs(A(k:end,k))); 
     r = n-(n-k+1)+r;  

     A([k r],:) = A([r k],:); 
     P([k r],:) = P([r k],:); 
     L([k r],:) = L([r k],:); 

     % from the pivot down divide by the pivot 
     L(k:n,k) = A(k:n,k)/A(k,k); 

     U(k,1:n) = A(k,1:n); 
     A(k+1:n,1:n) = A(k+1:n,1:n) - L(k+1:n,k)*A(k,1:n); 

    end 
    U(:,end) = A(:,end); 

end