LU分解
1.LU分解的基础知识
矩阵的LU分解又称为矩阵的三角分解，即将一个矩阵分解为一个下三角矩阵L和一个上三角矩阵U，即
A
=
L
U
A=LU
A=LU,其在方程组的求解和求矩阵的逆有许多应用。
LU分解的求解命令是lu，基本使用格式如下：
| 调用格式 |
说明 |
| [L,U]=lu(A) |
将一个矩阵分解为一个下三角L和一个上三角U |
| [L,U,P]=lu(A) |
对A进行分解，其中L为单位下三角矩阵，U为上三角矩阵，P为置换矩阵，满足LU=PA |
例1,对矩阵
A
=
[
1
2
3
4
5
6
7
8
2
3
4
1
7
8
5
6
]
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=⎣⎢⎢⎡1527263837454816⎦⎥⎥⎤进行分解。
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进行分解，得到
L
U
=
P
A
LU=PA
LU=PA,然后利用
L
y
=
P
b
Ly=Pb
Ly=Pb,最后再解
U
x
=
y
Ux=y
Ux=y,从而求解方程组的解。
例2：利用LU分解法求解方程组
{
x
1
+
x
2
−
3
x
3
−
x
4
=
1
3
x
1
−
x
2
−
3
x
3
+
4
x
4
=
4
x
1
+
5
x
2
−
9
x
3
−
8
x
4
=
0
\left\{\begin{array}{l}x_{1}+x_{2}-3 x_{3}-x_{4}=1 \\ 3 x_{1}-x_{2}-3 x_{3}+4 x_{4}=4 \\ x_{1}+5 x_{2}-9 x_{3}-8 x_{4}=0\end{array}\right.
⎩⎨⎧x1+x2−3x3−x4=13x1−x2−3x3+4x4=4x1+5x2−9x3−8x4=0 的解。
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
