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

合集下载
  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

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

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

相关文档
最新文档