《数值分析(06)矩阵分解法》由会员分享，可在线阅读，更多相关《数值分析(06)矩阵分解法（69页珍藏版）》请在金锄头文库上搜索。
1、,顺序高斯消元的基本思想：将矩阵A的下三角部分消为零，即,第三节 矩阵的三角分解法 一、矩阵三角分解的基本定理 顺序高斯消元与LU分解的等价性,进行到第k步消元时,其中,消元过程全部完成后，原来的二维数组中存放的元素实际上是一个新的矩阵，记为,function A=lud(A) %功能：对方阵A作三角分解A=LU,其中， % L为单位下三角阵，U为上三角阵， %输入：方阵A。 %输出：紧凑存储A=LU. %注意：当A的主元=0时退出Matlab.,LU分解的MATLAB程序,n,n=size(A); % 确定A的维数,for k=1:n-1for i=k+1:n if A(k,k) =0 qu
2、it; endA(i,k) =A(i,k)/ A(k,k);A(i,k+1:n)= A(i,k+1:n)- A(i,k) *A(k,k+1:n);end end,function A=lud(A) %功能：对方阵A作三角分解A=LU,其中， % L为单位下三角阵，U为上三角阵， %输入：方阵A。 %输出：紧凑存储A=LU.,LU分解的MATLAB程序,n,n=size(A); % 确定A的维数,for k=1:n-1if A(k,k) =0 break; endfor i=k+1:n A(i,k) =A(i,k)/ A(k,k);A(i,k+1:n)= A(i,k+1:n)- A(i,k) *
3、A(k,k+1:n);end end,function A=lud(A) %功能：对方阵A作三角分解A=LU,其中， % L为单位下三角阵，U为上三角阵， %输入：方阵A。 %输出：紧凑存储A=LU. n,n=size(A); % 确定A的维数 for k=1:n-1if A(k,k) =0 break; endA(k+1:n,k) =A(k+1:n,k)/ A(k,k);A(k+1:n ,k+1:n)= A(k+1:n ,k+1:n)- A(k+1:n ,k) *A(k,k+1:n);end,LU分解的MATLAB程序,列主元LU分解,Lupd.m %功能：对方阵A作列主元三角分解PA=LU
4、,其中， % L为单位下三角阵，U为上三角阵，排列阵P % 用向量p表示。 %输入：方阵A。 %输出：紧凑存储LU=LU，以及p。 %注意：当A奇异时退出Matlab. function LU,p=lupd(A) %初始化 n=length(A); p=1:n; LU=A;,%分解过程 for k=1:n%搜索列主元iks,i=max(abs(LU(k:n,k); ik=i+k-1;,算法3-3 列主元三角分解：PA=LU,%判断矩阵的奇异性 if s=0 quit; end %行交换if ik=km=p(k); p(k)=p(ik); p(ik)=m;lk=LU(k,:); LU(k,:)=
5、LU(ik,:); LU(ik,:)=lk;end,%用消元法计算LU=LUif k=n break; endLU(k+1:n,k)=LU(k+1:n,k)/LU(k,k);LU(k+1:n,k+1:n)=LU(k+1:n,k+1:n)-LU(k+1:n,k)*LU(k,k+1:n); end,A=0 1 -1 2;-1 0 0 1;2 -1 0 1;1 3 0 1,L,U,P=lu(A),L,U=lu(A),MATLAB实现,lupdsv.m %功能：调用列主元三角分解函数LU,p=lupd(A) % 求解线性方程组Ax=b。 %解法：PA=LU, Ax=bPAx=Pb % LUx=Pb,
6、y=Ux % Ly=f=Pb, f(i)=b(p(i) %输入：方阵A，右端项b（行或列向量均可） %输出：解x（行向量）,function x=lupdsv(A,b) n=length(b); LU,p=lupd(A); y(1)=b(p(1); for i=2:ny(i)=b(p(i)-LU(i,1:i-1)*y(1:i-1); end x(n)=y(n)/LU(n,n); for i=(n-1):-1:1x(i)=(y(i)-LU(i,i+1:n)*x(i+1:n)/LU(i,i); end,全主元LU分解,Lupqd.m %功能：对方阵A作全主元三角分解PAQT=LU,其中， % L为
7、单位下三角阵，U为上三角阵，排列阵P % 和Q分别用向量p,q表示。 %输入：方阵A。 %输出：紧凑存储LU=LU，以及p和q。 %注意：当A奇异时退出Matlab. function LU,p,q=lupqd(A) %初始化 n=length(A); p=1:n; q=p; LU=A;,%分解过程 for k=1:n%搜索全主元(ik,jk) xk,I=max(abs(LU(k:n,k:n); %列最大值及所在行 s,j=max(xk);ik=I(j)+k-1;jk=j+k-1; %判断矩阵的奇异性 if s=0 quit; end %行交换和列交换if ik=km=p(k); p(k)=p
8、(ik); p(ik)=m;lk=LU(k,:); LU(k,:)=LU(ik,:); LU(ik,:)=lk;end,if jk=km=q(k); q(k)=q(jk); q(jk)=m;ck=LU(:,k); LU(:,k)=LU(:,jk); LU(:,jk)=ck; end %用消元法计算LU=LUif k=n break; endLU(k+1:n,k)=LU(k+1:n,k)/LU(k,k);LU(k+1:n,k+1:n)=LU(k+1:n,k+1:n)-LU(k+1:n,k)*LU(k,k+1:n); end,lupqdsv.m %功能：调用全主元三角分解函数LU,p,q=lupq
9、d(A) % 求解线性方程组Ax=b。 %解法：PAQ-1=LU, Ax=b(PAQ-1)(Qx)=Pb % LU(Qx)=Pb, z=Qx, y=Uz % Ly=f=Pb, f(i)=b(p(i) % Uz=y, z=Qx , x(q(i)=z(i). %输入：方阵A，右端项b（行或列向量均可） %输出：解x（行向量）,function x=lupqdsv(A,b) n=length(b); LU,p,q=lupqd(A); y(1)=b(p(1); for i=2:ny(i)=b(p(i)-LU(i,1:i-1)*y(1:i-1); end z(n)=y(n)/LU(n,n);x(q(n)
10、=z(n); for i=(n-1):-1:1z(i)=(y(i)-LU(i,i+1:n)*z(i+1:n)/LU(i,i);x(q(i)=z(i); end,1.对称正定阵的Cholesky分解法,(1),三、矩阵的直接三角分解法,Cholesky分解,(k),若已求出了L的前k-1列元素，则为求第k列元素， 先用L的第k行乘 LT的第k列,再用L的第k行乘 LT的第j列( j=k+1 ,k+2,n )有,Cholesky分解公式,不完全Cholesky分解,定义 若n 阶矩阵A=(aij)的元素满足:对于1p,q a=-1*ones(1,4); b=4*ones(1,4); c=-1*on
11、es(1,4);d=3 2 2 3; x=zg(a,b,c,d) x = 1 1 1 1,周期三对角方程组的一般形式,基本思想：利用谢尔曼-莫里森公式（Sherman-Morrison)将方程组化为三对角方程组求解。,3、周期三对角方程组,谢尔曼-莫里森公式（Sherman-Morrison),如何选取U,V,习题三 P120-1 , 9,，使,、,寻找,。,a,a,j,i,k,ij,n,j,i,k,k,j,i,k,k,k,k,|,|,max,|,|,),(,),(,=,编程：,习题,，使,、,寻找,。,a,a,j,i,k,ij,n,j,i,k,k,j,i,k,k,k,k,|,|,max,|,|,),(,),(,=,编程：,
