当前位置:文档之家› 矩阵的LU分解(自编MATLAB)实验报告

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

LU 分解原理

定理:设ACnn,如果 A 的顺序主子式

𝐴11≠0, |𝑎11𝑎12𝑎21𝑎22|≠0,…,|𝑎11𝑎12𝑎21𝑎22…𝑎12…𝑎22⋮⋮𝑎𝑛−11𝑎𝑛−12⋮⋯𝑎𝑛−1𝑛−1|≠0

则存在唯一的主对角线上元素全为 1 的下三角矩阵L与唯一的上三角矩阵 U,使得

A=LU.

证明:对矩阵A的阶数使用数学归纳法.

显然,当 n=1 时,𝐴11=1 ∙𝐴11 就是唯一的分解式。现假定对

n-1 阶矩阵,定理的结论成立。对 A 进行分块

A=(𝐴𝐴−𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴)

其中𝐴𝐴,𝐴𝐴∈𝐴𝐴−𝐴.由于 n-1 阶矩阵 𝐴𝐴−𝐴的 k 阶顺序主子式就是 A 的 k 阶主子式(k=1,2,…,n-2),故它们都不为零.从而由归纳法假设,𝐴𝐴−𝐴 有唯一的 LU 分解

𝐴𝐴−𝐴=𝐴𝐴−𝐴𝐴𝐴−𝐴

其中𝐴𝐴−𝐴的主对角线上的元素都1.由于

|𝐴𝐴−𝐴|=|𝐴11𝐴12𝐴21𝐴22…𝐴12…𝐴22⋮⋮𝐴𝐴−11𝐴𝐴−12⋮⋯𝐴𝐴−1𝐴−1|=|𝐴𝐴−𝐴𝐴𝐴−𝐴|≠0

所以𝐴𝐴−𝐴及𝐴𝐴−𝐴是n-1阶可逆矩阵

先假设已有 A=LU,其中

L=(𝐴𝐴−𝐴0𝐴𝐴1), U= (𝐴𝐴−𝐴𝐴𝐴𝐴𝐴𝐴𝐴)

𝐴,𝐴∈𝐴𝐴−𝐴是待定向量。作乘积

𝐴𝐴 = (𝐴𝐴−𝐴𝐴𝐴−𝐴𝐴𝐴−𝐴𝐴𝐴𝐴𝐴𝐴−𝐴𝐴𝐴𝐴+𝐴𝐴𝐴) =(𝐴𝐴−𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴)=A

则𝐴,𝐴必须满足

𝐴𝐴−𝐴𝐴=𝐴𝐴,𝐴𝐴𝐴𝐴−𝐴=𝐴𝐴𝐴,𝐴𝐴𝐴+𝐴𝐴𝐴=𝐴𝐴𝐴

注意到𝐴𝐴−𝐴及𝐴𝐴−𝐴都是n-1阶可逆矩阵,则由上式可惟一确定

𝐴=𝐴𝐴−𝐴−𝐴𝐴𝐴,𝐴𝐴=𝐴𝐴𝐴𝐴𝐴−𝐴−𝐴, 𝐴𝐴𝐴=𝐴𝐴𝐴−𝐴𝐴𝐴

这就证明了 A 的 LU 分解的存在性和唯一性.

LU分解算法

当 n 阶矩阵满足定理的条件时,可以用初等变换的方法求出 L

和 U.

因为当 A=LU 时,由于 L 可逆,故必存在可逆矩阵 P 使得

𝐴𝐴=𝐴

即 PA=PLU=U.也就是说,可以先对 A 施行行的初等变换得出上三角矩阵U,而矩阵 P可以通过对单位矩阵I进行相同的行初等变换得出,即

P(A,I) =(PA,PI) =(U,P)

于是𝐴=𝐴−𝐴𝐴,为保持P为下三角矩阵(从而𝐴−𝐴也是下三角矩阵),在进行行初等变换时,不能进行行的对换,上行的倍数应加到下行的对应元.

LU分解用于解方程组

矩阵的三角分解在求解线性方程组时十分方便.如对线性方程组

𝐴𝐴=𝐴,设𝐴=𝐴𝐴.我们先求解方程组 𝐴𝐴=𝐴. 由于𝐴是下三角矩阵,则解向量𝐴可以通过依次求出其分量 𝐴1,𝐴2,⋯𝐴𝐴而求出,在求解方程组𝐴𝐴=𝐴.解向量𝐴可以通过该方程组依次求出分量𝐴𝐴,⋯,𝐴2,𝐴1而快速得出.于是由两个方程组𝐴𝐴=𝐴,𝐴𝐴=𝐴的求解而给出

𝐴𝐴𝐴=𝐴𝐴=𝐴= 𝐴𝐴的解.

程序流程图

输入矩阵A判断A是否为n×n矩阵?计算A的n-1阶顺序主子式是否为0?输出:“无法进行LU分解”设定 n 阶单位矩阵L和全零矩阵 U输出L,U是是否否计算U 的第一行

U1j =A1j j=1,2,…,n计算L 的第一列U 的第 r

行(逐行算出)L 的第 r

列(逐列算出)

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))

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

实际数据计算

已知矩阵A=(211410−221),𝐴=(121),先对A进性LU分解,并求解方程 A 𝐴=𝐴的解.

(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

0

2

x =

数据分析

调用MATLAB固有的LU分解函数,以及解方程组相关函数对以上数据进行计算,运行结果如下:

>> A=[2 1 1;4 1 0;-2 2 1];

>> b=[1 2 1]';

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

L =

0 0

0

U =

0

0

0 0

>> X=inv(A)*b

X =

经比对结果相同,可见以上程序可行。

相关主题