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

《矩阵的LU分解(自编MATLAB)实验报告.pdf》由会员分享，可在线阅读，更多相关《矩阵的LU分解(自编MATLAB)实验报告.pdf（7页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、1 1 矩阵的矩阵的 LULU 分解分解1.1LU1.1LU 分解原理分解原理定理：设 A AC11110, |21，如果 A A 的顺序主子式11121221222212| 022| 0,，|111211nn则存在唯一的主对角线上元素全为 1 的下三角矩阵 L L 与唯一的上三角矩阵 U U，使得A A=LULU. .证明：证明：对矩阵 A A 的阶数使用数学归纳法.显然， 当 n=1 时， 11=1 11就是唯一的分解式。 现假定对 n-1阶矩阵，定理的结论成立。对 A A 进行分块A A=()其中, .由于 n-1 阶矩阵 的 k 阶顺序主子式就是 A A 的 k 阶主子式k=1,2，n

2、-2 ，故它们都不为零.从而由归纳法假设，有唯一的 LULU 分解= 其中的主对角线上的元素都 1.由于111212212222|=|=| 0111211所以及是 n-1 阶可逆矩阵先假设已有 A A=LULU，其中L L=( = ( 则，必须满足0)，U=U=(1)=A A， 是待定向量。作乘积) =(+ = ，= ，+ = 注意到及都都是 n-1 阶可逆矩阵，则由上式可惟一确定 = ， = , , = 这就证明了 A A 的 LULU 分解的存在性和唯一性.1.2LU1.2LU 分解算法分解算法当 n 阶矩阵满足定理的条件时，可以用初等变换的方法求出 L L和 U U.因为当 A=LUA=

3、LU 时，由于 L L 可逆，故必存在可逆矩阵 P P 使得 = 即 PA=PLU=UPA=PLU=U.也就是说，可以先对 A A 施行行的初等变换得出上三角矩阵U U，而矩阵 P P可以通过对单位矩阵I I进行相同的行初等变换得出，即P(A,I)P(A,I)= (PA,PI)(PA,PI)= (U,P)(U,P)于是 = ,为保持P P为下三角矩阵从而也是下三角矩阵 ，在进行行初等变换时，不能进行行的对换，上行的倍数应加到下行的对应元.1.3LU1.3LU 分解用于解方程组分解用于解方程组矩阵的三角分解在求解线性方程组时十分方便.如对线性方程组 = ,设 = .我们先求解方程组 = . 由于

4、是下三角矩阵，则解向量可以通过依次求出其分量1,2,而求出，在求解方程组 = .解向量可以通过该方程组依次求出分量 ,2,1而快速得出.于是由两个方程组 = ， = 的求解而给出 = = = 的解.程序流程图程序流程图输入矩阵A判断A是否为nn矩阵？是否计算A的n-1阶顺序主子式是否为0？否是设定 n 阶单位矩阵L和全零矩阵 U计算U 的第一行U1j=A1j j=1，2，n计算L 的第一列U 的第 r行（逐行算出）L 的第 r列（逐列算出）输出L,U输出：“无法进行LU分解”1.5MATLAB1.5MATLAB 程序程序function f=LU_decom(A)m,n=size(A)if m

5、=nfprintf(Error:m and n must be equal!m=%d,n=%dn,m,n)endfor i=1:n-1if (det(A(1:i,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)=

6、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 =100210-1-31U =2110-1-200-42解方程组，程序及结果如下%-用LU分解解线性方程组-y=zeros(n,1);y(1)=b(1);for i=2:ny(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:1x(i)=(y(i)-sum(U(i,i+1:n).*x(i+1)/U(i,i);endx=x运行结果如下：y =102x =数据分析数据分析调用MATLAB固有的LU分解函数，以及解方程组相关函数对以上数据进行计算，运行结果如下： A=2 1 1;4 1 0;-2 2 1; b=1 2 1; L,U=lu(A)L =1.000000-0.50001.00000U =4.00001.00000 X=inv(A)*bX =经比对结果相同，可见以上程序可行。