《矩阵的LU分解(自编MATLAB)实验报告(总7页)》由会员分享，可在线阅读，更多相关《矩阵的LU分解(自编MATLAB)实验报告(总7页)（7页珍藏版）》请在金锄头文库上搜索。
1、1矩阵的LU分解1.1 LU 分解原理定理：设ACnn，如果 A 的顺序主子式A110, a11a12a21a220,，a11a12a21a22a12a22an-11an-12an-1n-10则存在唯一的主对角线上元素全为 1 的下三角矩阵L与唯一的上三角矩阵 U，使得A=LU.证明：对矩阵A的阶数使用数学归纳法.显然，当 n=1 时，A11=1 A11 就是唯一的分解式。现假定对 n-1 阶矩阵，定理的结论成立。对 A 进行分块A=An-112Tnn其中1,2Cn-1.由于 n-1 阶矩阵 An-1的 k 阶顺序主子式就是 A 的 k 阶主子式（k=1,2，n-2），故它们都不为零.从而由归
2、纳法假设，An-1 有唯一的 LU 分解An-1=Ln-1Un-1其中Ln-1的主对角线上的元素都1.由于An-1=a11a12a21a22a12a22an-11an-12an-1n-1=Ln-1Un-10所以Ln-1及Un-1是n-1阶可逆矩阵先假设已有 A=LU，其中 L=Ln-10T1， U= Un-1Tbnn，Cn-1是待定向量。作乘积LU = Ln-1Un-1Ln-1TUn-1bnn+T =An-112Tnn=A则，必须满足 Ln-1=1，TUn-1=2T，bnn+T=nn注意到Ln-1及Un-1都是n-1阶可逆矩阵，则由上式可惟一确定=Ln-1-11，T=2TUn-1-1, bnn
3、=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分解用于解方程组矩阵的三角分解在求解线性方程组
4、时十分方便.如对线性方程组Ax=b,设A=LU.我们先求解方程组 Ly=b. 由于L是下三角矩阵，则解向量y可以通过依次求出其分量 y1,y2,yn而求出，在求解方程组Ux=y.解向量x可以通过该方程组依次求出分量xn, ,x2,x1而快速得出.于是由两个方程组Ux=y，Ly=b的求解而给出LUx=Ly=b= Ax的解.1.4程序流程图1.5 MATLAB程序function f=LU_decom(A) m,n=size(A)if m=nfprintf(Error:m and n must be equal!m=%d,n=%dn,m,n) endfor i=1:n-1if (det(A(1:i
5、,1:i)=0)fprintf(Error:det A(%d,%d)=0!n,i,i)flag=failurereturn;elseflag=ok;endendL=eye(n);U=zeros(n);for i=1:nU(1,i)=A(1,i);endfor r=2:nL(r,1)=A(r,1)/U(1,1);endfor i=2:nfor j=i:nz=0;for r=1:i-1z=z+L(i,r)*U(r,j);endU(i,j)=A(i,j)-z;endif abs(U(i,i) A=2 1 1;4 1 0;-2 2 1;LU_decom(A)m = 3n = 3L = 1 0 0 2
6、1 0 -1 -3 1U = 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);endyx(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);endx=x运行结果如下：y = 1 0 2x = -0.5000 1.0000 -0.50001.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 0U = 4.0000 1.0000 0 0 2.5000 1.0000 0 0 0.8000 X=inv(A)*bX = 0.2500 1.0000 -0.5000经比对结果相同，可见以上程序可行。
