高质量学习资源免费获取，专注但不限于【Linux】【C/C++/Qt】【FPGA】【数据结构与算法】, 根据多年技术经验纯【原创】，纯【干货】，分享【技术心得】，力求【授人以鱼，更授人以渔】。

　　高斯消元法求解线性方程，包括把增广矩阵转换为三角矩阵形式的过程，消去阶段工作

步骤是把矩阵A分解成为下三角L和上三角U的乘积。这种计算L,U的过程称为LU分解法。

lu实现对矩阵的分解。

[L,U] = lu(A)

%%将矩阵A分解的上三角矩阵保存在U当中，将一个“心理学上的”下三角矩阵（例如一个下三角矩阵和置换矩阵的乘积）保存在L中，满足A=L*U，注意A不必须是方阵。

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

%%返回三个矩阵，下三角矩阵L、上三角矩阵U和一个置换矩阵P，满足P*A=L*U。

[L,U,p] = lu(A,'vector')

%%将使用向量而不是矩阵的方式返回置换信息，p是一个行向量，满足A(p,:) = L*U，如果使用[L,U,P] = lu(A,'matrix')则返回的P是向量，这也是默认的形式。

[L1,U1]=lu(x)
[L2,U2,P]=lu(x)
[L3,U3,P,Q] = lu(X)
MATLAB中[L1,U1]=lu(x)的结果：L是下三角的置换矩阵即L1=p*L2，U是上三角阵。Matlab中LU分解采用高斯消元法，结果是不唯一的，只要[L1,U1]=lu(x)满足L1*U1=x, [L2,U2,P]=lu(x)满足L2*U2=p*x，[L3,U3,P,Q] = lu(X)满足 L3*U3= P*X*Q就行了

例方程组:

x1+0.2*(x2)+0.5*(x3)=1

0.2*(x1)+0.4*(x2)+x3=2

0.5*(x1)+0.1*(x2)+0.6*(x3)=1.5

根据LU分解有LUx=b，记Y=L^(-1)b，则x=U^(-1)Y为方程组Ax=b解。

>> a=[1 0.2 0.5;0.2 0.4 1; 0.5 0.1 0.6];
>> b=[1 2 1.5];
>> [L,U]=lu(a)

L =

    1.0000         0         0
    0.2000    1.0000         0
    0.5000         0    1.0000

U =

    1.0000    0.2000    0.5000
         0    0.3600    0.9000
         0         0    0.3500

>> y=inv(L)*b'  

y =

    1.0000
    1.8000
    1.0000

或者 

>> y=L^(-1)*b'

y =

    1.0000
    1.8000
    1.0000

>> x=U^(-1)*y

x =

         0
   -2.1429
    2.8571

或

>> x=inv(U)*y

x =

         0
   -2.1429
    2.8571

下面就是lu函数的厂商代码，相信这段程序，对大家编写线性方程组直接解法的诸多算法都有启发作用。

 1 function [L,U,x]=lux(A,b)
 2 %LU 分解法解线性方程组（列主元LU分解）
 3 [n,n]=size(A);
 4 p=eye(n);%p记录了选择主元时候所进行的行变换
 5 for k=1:n-1
 6      [r,m]=max(abs(A(k:n,k)));  %选列主元
 7      m=m+k-1;
 8      if(A(m,k)~=0)
 9          if(m~=k)
10              A([k m],:)=A([m k],:);
11              p([k m])=p([m k]);
12         end
13         for i=k+1:n    
14              A(i,k)=A(i,k)/A(k,k);
15              j=k+1:n;              A(i,j)=A(i,j)-A(i,k)*A(k,j);
16         end
17     end
18 end
19 L=tril(A,-1)+eye(n,n);
20 U=triu(A);
21 %解下三角矩阵 Ly=b
22  newb=p*b;
23  y=zeros(n,1);
24 for k=1:n
25     j=1:k-1;
26     y(k)=(newb(k)-L(k,j)*y(j))/L(k,k);
27 end 
28 %解上三角方程组  Ux=y
29 x=zeros(n,1);
30 for k=n:-1:1
31      j=k+1:n;
32      x(k)=(y(k)-U(k,j)*x(j))/U(k,k);
33 end