線性方程組的分解法——LU分解法

　　1.代碼

%%LU分解法
function LUDM = LU_Decomposition_method(A,b)
global n;global B;global U;global L;global M;
[n,n] = size(A);
B = [A,b];
R_A = rank(A);R_B = rank(B);
if R_A ~= R_B
    disp('方程無解');
elseif (R_A == R_B) && (R_A == n)
    disp('此方程有唯一解');
    M = LU_decomposition(A);
    L = M(:,:,1);U = M(:,:,2);
    matrix1 = [L b];
    Y = Lower_trig_iterative_solution(matrix1);
    matrix2 = [U Y];
    X = Upper_trig_iterative_solution(matrix2);
    disp('LU分解中L=');
    L
    disp('LU分解中U=');
    U
else
    disp('方程有無窮多組解');
end
disp('解向量為:');
LUDM = X;

%%矩陣的LU分解
    function LUD = LU_decomposition(A)
        [n,n] = size(A);
        M = Elementary_transformation_of_the_lower_triangle(A);
        L = M(:,:,n);U=A;
        for i = 1:1:n-1
            U = M(:,:,i)*U;
        end
        LUD(:,:,1) = L;
        LUD(:,:,2) = U;
    end
%%下三角初等變換
    function ETLT = Elementary_transformation_of_the_lower_triangle(A)
        [n,n] = size(A);
        L = zeros(n,1,n);
        for i = 1:1:n
            for j = 1:1:n
                for k = 1:1:n
                    if j == k
                        L(j,k,i) = 1;
                    end
                end
            end
        end
        for i = 1:1:n-1
            for j = 1:1:n
                for k = 1:1:n
                    if j > k
                        if i == k
                            L(j,k,i) = -A(j,k)/A(k,k);
                        end
                        L(i+1:n,i,n) = -L(i+1:n,i,i);
                    end
                end
            end
            A = L(:,:,i)*A;
        end
        ETLT = L;
    end
%%下三角迭代法
    function LTIS = Lower_trig_iterative_solution(M)
        [m,n] = size(M);
        B  =M(:,1:n-1);ba = M(:,n);
        y = zeros(1,m);
        y(1) = ba(1);
        for i = 2:1:m
            sum = 0;
            for j = 1:1:i-1
                sum = sum+B(i,j)*y(j);
            end
            y(i) = ba(i)-sum;
        end
        LTIS = y';
    end
%%上三角迭代法
    function UTIS = Upper_trig_iterative_solution(M)
        [m,n] = size(M);
        B = M(:,1:n-1);ba = M(:,n);
        x = zeros(1,m);
        x(m) =ba(m)/B(m,m);
        for i = m-1:-1:1
            sum = 0;
            for j = i+1:1:m
                sum = sum+B(i,j)*x(j);
            end
            x(i) = (ba(i)-sum)/B(i,i);
        end
        UTIS = x';
    end
end

　　2.例子

clear all
clc
M = rand(9)
b = reshape(rand(3),9,1)
 
S = LU_Decomposition_method(M,b)

M\b

　　結果

M =
  列 1 至 7
    0.5944    0.4709    0.4076    0.4235    0.5181    0.0680    0.6022
    0.0225    0.6959    0.8200    0.0908    0.9436    0.2548    0.3868
    0.4253    0.6999    0.7184    0.2665    0.6377    0.2240    0.9160
    0.3127    0.6385    0.9686    0.1537    0.9577    0.6678    0.0012
    0.1615    0.0336    0.5313    0.2810    0.2407    0.8444    0.4624
    0.1788    0.0688    0.3251    0.4401    0.6761    0.3445    0.4243
    0.4229    0.3196    0.1056    0.5271    0.2891    0.7805    0.4609
    0.0942    0.5309    0.6110    0.4574    0.6718    0.6753    0.7702
    0.5985    0.6544    0.7788    0.8754    0.6951    0.0067    0.3225
  列 8 至 9
    0.7847    0.1917
    0.4714    0.7384
    0.0358    0.2428
    0.1759    0.9174
    0.7218    0.2691
    0.4735    0.7655
    0.1527    0.1887
    0.3411    0.2875
    0.6074    0.0911
b =
    0.5762
    0.6834
    0.5466
    0.4257
    0.6444
    0.6476
    0.6790
    0.6358
    0.9452
此方程有唯一解
LU分解中L=
L =
  列 1 至 7
    1.0000         0         0         0         0         0         0
    0.0379    1.0000         0         0         0         0         0
    0.7155    0.5352    1.0000         0         0         0         0
    0.5261    0.5762  -74.4491    1.0000         0         0         0
    0.2717   -0.1391 -136.4397    1.7669    1.0000         0         0
    0.3008   -0.1074  -74.0359    0.9200    0.6765    1.0000         0
    0.7115   -0.0228   42.5434   -0.5996    0.3838 -141.0829    1.0000
    0.1585    0.6728   -1.3001   -0.0414    0.8852  -70.1396    0.4925
    1.0070    0.2658  -39.5864    0.4476    1.3552   49.3425   -0.3788
  列 8 至 9
         0         0
         0         0
         0         0
         0         0
         0         0
         0         0
         0         0
    1.0000         0
    5.1107    1.0000
LU分解中U=
U =
  列 1 至 7
    0.5944    0.4709    0.4076    0.4235    0.5181    0.0680    0.6022
         0    0.6781    0.8045    0.0748    0.9240    0.2522    0.3640
         0         0   -0.0039   -0.0765   -0.2275    0.0404    0.2903
         0         0         0   -5.8101  -16.7848    3.4944   21.0900
   -0.0000         0         0         0   -1.1550    0.1988    2.6992
    0.0000         0         0         0         0   -0.0074    0.5483
    0.0000   -0.0000         0         0         0         0   76.6535
    0.0000    0.0000         0   -0.0000         0         0         0
   -0.0000   -0.0000         0    0.0000         0         0         0
  列 8 至 9
    0.7847    0.1917
    0.4416    0.7312
   -0.7621   -0.2857
  -57.2283  -20.8735
   -2.2924   -1.7782
   -1.9343    0.0429
 -274.3037    6.4447
   -1.9999   -0.0598
         0    0.7768
解向量為:
S =
   -0.9496
    2.2130
    0.5483
    1.9595
   -3.8859
   -0.4632
    0.4453
    0.3978
    2.6573
ans =
   -0.9496
    2.2130
    0.5483
    1.9595
   -3.8859
   -0.4632
    0.4453
    0.3978
    2.6573
>>