MATLAB之LU分解法（十）

LU分解

1.LU分解的基础知识

矩阵的LU分解又称为矩阵的三角分解，即将一个矩阵分解为一个下三角矩阵L和一个上三角矩阵U，即$A=LU$,其在方程组的求解和求矩阵的逆有许多应用。

LU分解的求解命令是lu，基本使用格式如下：

调用格式	说明
[L,U]=lu(A)	将一个矩阵分解为一个下三角L和一个上三角U
[L,U,P]=lu(A)	对A进行分解，其中L为单位下三角矩阵，U为上三角矩阵，P为置换矩阵，满足LU=PA

例1,对矩阵$A=\left[\begin{array}{llll}1& 2& 3& 4\\ 5& 6& 7& 8\\ 2& 3& 4& 1\\ 7& 8& 5& 6\end{array}\right]$进行分解。

A=[1,2,3,4;5,6,7,8;2,3,4,1;7,8,5,6];
[L1,U1]=lu(A)

L1 =

​    0.1429    1.0000         0         0
​    0.7143    0.3333    1.0000         0
​    0.2857    0.8333    0.2500    1.0000
​    1.0000         0         0         0


U1 =

​    7.0000    8.0000    5.0000    6.0000
​         0    0.8571    2.2857    3.1429
​         0         0    2.6667    2.6667
​         0         0         0   -4.0000

[L2,U2,P]=lu(A)

L2 =

​    1.0000         0         0         0
​    0.1429    1.0000         0         0
​    0.7143    0.3333    1.0000         0
​    0.2857    0.8333    0.2500    1.0000


U2 =

​    7.0000    8.0000    5.0000    6.0000
​         0    0.8571    2.2857    3.1429
​         0         0    2.6667    2.6667
​         0         0         0   -4.0000


P =

 0     0     0     1
 1     0     0     0
 0     1     0     0
 0     0     1     0

从上面两者的计算结果比较分析可以得出：使用第一种方法得到的上三角矩阵并不是严格的上三角矩阵，这样对于以后的分析和计算显然是不利的，所以在工程计算或者其他计算中都采用第二种方法。

2.LU分解法求解方程组

LU分解法求解的基本思路就是：将系数矩阵A进行分解，得到$LU=PA$,然后利用$Ly=Pb$,最后再解$Ux=y$,从而求解方程组的解。

例2：利用LU分解法求解方程组$\{\begin{array}{l}{x}_1+x_2-3x_3-x_4=1\\ 3x_1-x_2-3x_3+4x_4=4\\ x_1+5x_2-9x_3-8x_4=0\end{array}$ 的解。

A = [1,1,-3,-1;3,-1,-3,4;1,5,-9,-8];
b=[1;4;0]

b =

     1
     4
     0

[L,U,P]=lu(A)

L =

    1.0000         0         0
    0.3333    1.0000         0
    0.3333    0.2500    1.0000


U =

    3.0000   -1.0000   -3.0000    4.0000
         0    5.3333   -8.0000   -9.3333
         0         0         0    0.0000


P =

     0     1     0
     0     0     1
     1     0     0

y = L'*P*b

y =

    4.3333
    0.2500
    1.0000

>> x = U'*y

x =

   13.0000
   -3.0000
  -15.0000
   15.0000

利用函数实现如下：

function x= solvebyLU(A,b)
%此函数用于利用LU分解法解方程组
flag=isexist(A,b);    %判断方程组是否有解
if flag==0
    disp('方程组无解');
    x = [];
    return;
else
    r = rank(A);
    [m,n]=size(A);
    [L,U,P]=lu(A);
    b = P*b;
    

    %解Ly=b
    y(1) = b(1);
    if m>1
        for i=2:m
            y(i)=b(i)-L(i,1:i-1)*y(1:i-1)';
        end
    end
    y = y';
    
    %解Ux=y得方程组的解
    x0(r)=y(r)/U(r,r);
    if r>1
        for i=r-1:-1:1
            x0(i)=(y(i)-U(i,i+1:r)*x0(i+1:r)')/U(i,i);
        end
    end
    x0 = x0';
    
    if flag==1
        x=x0;
        return;
    else
        format rat;
        Z = null(A,'r');
        [mz,nz]=size(Z);
        x0(r+1:n)=0;
        for i=1:nz
            t=sym(char([107 48+i]));
            k(i)=t;
        end
        x=x0;
        for i=1:nz
            x = x+k(i)*Z(:,i);
        end
    end

end

运行结果：

 A = [1,1,-3,-1;3,-1,-3,4;1,5,-9,-8];

b = [1;4;0];
x = solvebyLU(A,b)

x =

 (3*k1)/2 - (3*k2)/4 + 5/4
 (3*k1)/2 + (7*k2)/4 - 1/4
                        k1
                        k2