5. 直接三角形分解法

MATLAB函数文件

function [l,u,x,y] = agui_lu(a,b)
% 求可逆矩阵的分解，l为下三角矩阵，u为上三角矩阵
n = length(a);
u = zeros(n,n);
l = eye(n,n);
u(1,:)=a(1,:);
l(2:n,1)=a(2:n,1)/u(1,1);
for k=2:n
    u(k,k:n)=a(k,k:n)-l(k,1:k-1)*u(1:k-1,k:n);
    l(k+1:n,k)=( a(k+1:n,k)-l(k+1:n,1:k-1)*u(1:k-1,k) )/u(k,k);
end
l
u
% 解Ly=b
y = zeros(n,1);
y(1)=b(1);
for k=2:n
    y(k)=b(k)-l(k,1:k-1)*y(1:k-1);
end
y
%  解Ux=y
x = zeros(n,1);
x(n)=y(n)/u(n,n);
for k=(n-1):-1:1
    x(k)=(y(k)-u(k,(k+1):n)*x((k+1):n))/u(k,k);
end
x
end

测试

>> a=[2 1 5;4 1 12;-2 -4 5]

a =

     2     1     5
     4     1    12
    -2    -4     5

>> b=[11;27;12]

b =

    11
    27
    12

>> agui_lu(a,b)

l =

     1     0     0
     2     1     0
    -1     3     1


u =

     2     1     5
     0    -1     2
     0     0     4


y =

    11
     5
     8


x =

     1
    -1
     2