矩阵的LU分解(自编MATLAB)实验报告

1矩阵的LU分解
1.1 LU 分解原理
定理：设A C n n，如果A的顺序主子式
A11≠0, |a11a12
a21a22|≠0,…，|
a11a12
a21a22
…a12
…a22
⋮⋮
a n−11a n−12
⋮
⋯a n−1n−1
|≠0
则存在唯一的主对角线上元素全为1 的下三角矩阵L与唯一的上三角矩阵U，使得
A=LU.
证明：对矩阵A的阶数使用数学归纳法.
显然，当n=1 时，A11=1 ∙A11就是唯一的分解式。

现假定对n-1 阶矩阵，定理的结论成立。

对A进行分块
A=(A n−1α1α2Tαnn
)
其中α1,α2∈C n−1.由于n-1 阶矩阵A n−1的k阶顺序主子式就是A的k阶主子式（k=1,2，…，n-2），故它们都不为零.从而由归纳法假设，A n−1有唯一的LU分解
A n−1=L n−1U n−1
其中L n−1的主对角线上的元素都1.由于
|A n−1|=|a11a12
a21a22
…a12
…a22⋮⋮
a n−11a n−12
⋮
⋯a n−1n−1
|=|L n−1U n−1|≠0
所以L n−1及U n−1是n-1阶可逆矩阵先假设已有A=LU，其中
L=(L n−10
γT1
)，U= (
U n−1β
γT b nn
)
β，γ∈C n−1是待定向量。

作乘积
LU =(L n−1U n−1L n−1β
γT U n−1b nn+γTβ
)=(
A n−1α1
α2Tαnn
)=A
则β，γ必须满足
L n−1β=α1，γT U n−1=α2T ，b nn +γT β=αnn
注意到L n−1及U n−1都是n-1阶可逆矩阵，则由上式可惟一确定
β=L n−1−1α1，γT =α2T U n−1−1, b nn =αnn −γT β
这就证明了 A 的 LU 分解的存在性和唯一性.
1.2 LU 分解算法
当 n 阶矩阵满足定理的条件时，可以用初等变换的方法求出 L 和 U .
因为当 A=LU 时，由于 L 可逆，故必存在可逆矩阵 P 使得
PL =I
即 PA=PLU=U .也就是说，可以先对 A 施行行的初等变换得出上三角矩阵U ，而矩阵 P 可以通过对单位矩阵I 进行相同的行初等变换得出，即
P(A,I) =(PA,PI) =(U,P)
于是A =P −1U ,为保持P 为下三角矩阵（从而P −1也是下三角矩阵），在进行行初等变换时，不能进行行的对换，上行的倍数应加到下行的对应元.
1.3 LU 分解用于解方程组
矩阵的三角分解在求解线性方程组时十分方便.如对线性方程组 Ax =b,设A =LU.我们先求解方程组 Ly =b. 由于L 是下三角矩阵，则解向量y 可以通过依次求出其分量 y 1,y 2,⋯y n 而求出，在求解方程组Ux =y.解向量x 可以通过该方程组依次求出分量x n ,⋯,x 2,x 1而快速得出.于是由两个方程组Ux =y ，Ly =b 的求解而给出
LUx =Ly =b = Ax 的解.
1.4程序流程图
1.5 MATLAB程序
function f=LU_decom(A)
[m,n]=size(A)
if m~=n
fprintf('Error:m and n must be equal!m=%d,n=%d\n',m,n) end
for i=1:n-1
if (det(A(1:i,1:i))==0)
fprintf('Error:det A(%d,%d)=0!\n',i,i)
flag='failure'
return;
else
flag='ok';
end
end
L=eye(n);
U=zeros(n);
for i=1:n
U(1,i)=A(1,i);
end
for r=2:n
L(r,1)=A(r,1)/U(1,1);
end
for i=2:n
for j=i:n
z=0;
for r=1:i-1
z=z+L(i,r)*U(r,j);
end
U(i,j)=A(i,j)-z;
end
if abs(U(i,i))<eps
flag='failure'
return;
end
for k=i+1:n
m=0;
for q=1:i-1
m=m+L(k,q)*U(q,i);
end
L(k,i)=(A(k,i)-m)/U(i,i); end end
L
U
1.6 实际数据计算
已知矩阵A=(211 410
−221)，b=(
1
2
1
)，先对A进性LU分解，
并求解方程A x=b的解.
（1）A的LU分解
在MATLAB命令行中输入A=[2 1 1;4 1 0;-2 2 1];并调用以上函数可得如下结果
>> A=[2 1 1;4 1 0;-2 2 1];LU_decom(A)
m =
3
n =
3
L =
1 0 0
2 1 0
-1 -3 1
U =
2 1 1
0 -1 -2
0 0 -4
（2）解方程组，程序及结果如下
%-----用LU分解解线性方程组------
y=zeros(n,1);
y(1)=b(1);
for i=2:n
y(i)=b(i)-sum(L(i,1:i-1)'.*y(1:i-1));
end
y
x(n)=y(n)/U(n,n);
for i=n-1:-1:1
x(i)=(y(i)-sum(U(i,i+1:n)'.*x(i+1)))/U(i,i);
end
x=x'
运行结果如下：
y =
1
2
x =
-0.5000
1.0000
-0.5000
1.7数据分析
调用MATLAB固有的LU分解函数，以及解方程组相关函数对以上数据进行计算，运行结果如下：
>> A=[2 1 1;4 1 0;-2 2 1];
>> b=[1 2 1]';
>> [L,U]=lu(A)
L =
0.5000 0.2000 1.0000
1.0000 0 0
-0.5000 1.0000 0
U =
4.0000 1.0000 0
0 2.5000 1.0000
0 0 0.8000
>> X=inv(A)*b
X =
0.2500
1.0000
-0.5000
经比对结果相同，可见以上程序可行。