中科院矩阵分析与应用大作业

精品文档中科院矩阵分析与应用大作业实现LU分解 QR分解 Householder reduction、Givens reduction Matlab 代码：function =juzhendazuoyeA=input(请输入一个矩阵A=);x=input(请输入序号 1 LU分解 2 Gram-Schmidt分解 3 Householder reduction 4 Givens reduction: ); if(x=1) %*LU分解*% disp(PA=LU)m=size(A,1); % m等于矩阵A的行数 n=size(A,2); % n等于矩阵A的列数if(m=n) % 判断矩阵A是不是方阵 % 如果矩阵A不是方阵那么就输出“error” U=A; % 把矩阵A赋值给矩阵UL=zeros(n); % 先将L设为单位阵P=eye(n); % 首先将交换矩阵P设为单位矩阵for j=1:n-1 for i=j+1:n if (U(j,j)=0) %判断主元元素是否不为0 L(i,j)=U(i,j)/U(j,j); U(i,:)=U(i,:)-U(j,:)*U(i,j)/U(j,j); % U(j,j)为主元元素 else a=j+1; % 令a等于j+1 while(U(a,j)=0)&(an) % 判断主元元素所对的下一行元素是不是0，a是否小于n a=a+1; % 寻找下一个元素 end temp=U(j,:); % 判断主元元素所在列（除主元元素外）第一个不为零的元素的所在行与主元元素所在行进行行交换 U(j,:)=U(a,:); % U两行交换位置 U(a,:)=temp ; m=L(j,:); L(j,:)=L(a,:); % L矩阵两行交换位置 L(a,:)=m; q=P(j,:); P(j,:)=P(a,:); % 交换矩阵的两行交换 P(a,:)=q; L(i,j)=U(i,j)/U(j,j); U(i,:)=U(i,:)-U(j,:)*U(i,j)/U(j,j); endendendfor k=1:n L(k,k)=1; % 把L矩阵的对角线赋值为1endL % 输出下三角矩阵L U % 输出上三角矩阵UP % 输出交换矩阵P A=inv(P)*L*Uelse disp(error)end end if(x=2) % 判断如果x=2，那么将执行schmid分解%*Gram-Schmidt正交分解*%disp(A=Q*R)Q=zeros(size(A,1),size(A,2); % 先把Q设为全零矩阵R=zeros(size(A,2); % R设置为全零矩阵a=A(:,1); % 把第一列赋值给aR(1,1)=norm(a); % 求第一列列向量的模值a=a/norm(a); % 求第一列列向量的单位向量Q(:,1)=a; % 把a赋值给Q的第一列for j=2:size(A,2) m=zeros(size(A,1),1); % 取A的第一列 for i=1:j-1 R(i,j)=Q(:,i)* A(:,j); % q的转置乘以A的第j列向量 m=m+R(i,j)*Q(:,i); % q的转置乘以A的列向量 end Q(:,j)=A(:,j)-m; % A的第j列减去q(i)和A(:,j)的内积和 R(j,j)=norm(Q(:,j); % 把Q的列向量的模值赋值给R(j,j) Q(:,j)=Q(:,j)/norm(Q(:,j); % 把Q的列向量的单位化end Q % 输出正交矩阵Q R % 输出上三角矩阵Rend if(x=3) % 判断如果x=3，那么将进行Householder reduction %*Householder reduction*%disp(P*A=T) R=zeros(size(A,1); % 把R设置为矩阵维数等于矩阵的行数的全零方阵R1=zeros(size(A,1); % 把R1设为矩阵维数等于矩阵的行数的全零方阵M=A; % 将A赋给M P=eye(size(M,1); % 先将P矩阵设为维数等于M的单位矩阵 for i=1:(size(M,1)-1) U=A; % 将A赋值给U U(1,1)=U(1,1)-norm(U(:,1); % 将U的第一列的第一行元素减去U的第一列列向量的模值 R=eye(size(U,1)-2*U(:,1)*U(:,1)/(U(:,1)* U(:,1); % I-2*U(:,1)*U(:,1)/(U(:,1)* U(:,1) A=R*A; % R乘以A赋值给A A=A(2:size(A,1),2:size(A,2); % 取A的子矩阵 if(size(R,1)size(M,1) % 判断矩阵R的行数是否小于矩阵M的行数，如果小于进行下步： S=eye(size(M,1)-size(R,1); % 将S设置为维数等于矩阵M的行数减去矩阵R的行数维的单位矩阵 V=zeros(size(M,1)-size(R,1),size(R,1); % 将V设置为矩阵行数等于M的行数减去R的行数，列数等于矩阵R的列数 F=zeros(size(R,1),size(M,1)-size(R,1); % 将F矩阵设置行数等于R的行数，列数等于矩阵M的行数减去矩阵R的行数 R1=S V;F R; % 将 S V F D 合成矩阵R1 else R1=R; % 如果不满矩阵R的行数小于矩阵M的行数，则把R赋值给R1 end P=R1*P; end P % 输出正交矩阵PT=P*M % 输出矩阵T,如果矩阵M的行数等于列数的话，T为上三角矩阵end if(x=4) % 判断x的值是否等于4，等于4则进行Givens reduction%*Givens reduction*%disp(P*A=R) U=A; % 将A赋值给Uw=size(A,1); % w等于矩阵A的行数 r=eye(w); % 将r设置为维数为w的单位矩阵for k=1:w-1 m=eye(size(A,1); % 将m设置为维数等于A的行数单位矩阵 for i=2:size(A,1) P=eye(size(A,1); a=0; % 将a是设置为0，方便求第一列前i个元素的平方和 for j=1:i u=sqrt(a); a=a+A(j,1)2; end s=sqrt(a); % 将第一列前i个元素的平方开根 P(1,1)=u/s; % 将u/s赋值给旋转矩阵P的第一行的第一列 P(i,i)=u/s; % 将u/s赋值给旋转矩阵P的第i行和第i列 P(i,1)=-A(i,1)/s; % 将 -A(i,1)赋值非P的第i行的第一列 P(1,i)=A(i,1)/s; % 将 A(i,i)赋值给P的第一行的第i列 m=P*m; % P乘以矩阵m并赋值给m endA=m*A; % 矩阵m*A赋值给AA=A(2:size(A,1),2:size(A,2); % 取A的子矩阵if(size(m,1)w) % 如果矩阵m的行数小于wc=eye(w-size(m,1); % 将c设置为维数等于w-矩阵m的行数的单位矩阵d=zeros(w-size(m,1),size(m,1);v=zero