MATLAB代码实现LU分解

LU分解

%一百行以下的小方阵
function [l,u] = DLU(a)
n = size(a);
n = n(1);
l=eye(n); 
u=zeros(n); 
for k=1:n 
    for j=k:n 
        b=0;
        for r=1:k-1
            b=b+l(k,r)*u(r,j);
        end
        u(k,j)=a(k,j)-b;
    end
    for i=k+1:n 
        b=0;
        for r=1:k-1 
            b=b+l(i,r)*u(r,k);
        end
        l(i,k)=(a(i,k)-b)/u(k,k);
    end
end

%一百行以上的大方阵
function [L,U] = DLU( A )
%DLU 返回对角元素全1的下三角矩阵L和上三角矩阵U，参数为一个n阶方阵
%
n=size(A);
n = n(1);
%开始计算LU矩阵
%生成n阶对角阵L
L = eye(n,n);
%生成n阶0矩阵U
U = zeros(n,n);
for i = 1:n-1
    %行规则减去内积
    A(i,i:end) = A(i,i:end) - A(i,1:i-1)*A(1:i-1,i:end);
    U(i,i:end) = A(i,i:end);
    %列规则减去内积除以对角元
    A(i+1:end,i) = (A(i+1:end,i) - A(i+1:end,1:i-1)*A(1:i-1,i))/A(i,i);
    L(i+1:end,i) = A(i+1:end,i);
end
%利用行规则计算最后一个元素
A(n,n) = A(n,n) - A(n,1:n-1)*A(1:n-1,n);
U(n,n) = A(n,n);
end

function FLAG = canLU( A )
%判断可不可以LU分解
FLAG = 0;
if ismatrix(A)==0
    disp('输入参数不是矩阵')
    return
end
[m,n] = size(A);
if m ~= n
    disp('矩阵不是方阵')
    return
end
for p = 1:n
    if rank(A(1:p,1:p)) ~= p
        disp('矩阵存在为零的顺序主子式，不能LU分解')
        return
    end
end
FLAG = 1;
end