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

中科院矩阵分析与应用大作业实现LU分解QR分解Householderreduction、GivensreductionMatlab代码：function[]=juzhendazuoyeA=input('请输入一个矩阵A=');x=input('请输入序号1LU分解2Gram-Schmidt分解3Householderreduction4Givensreduction:');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设为单位矩阵forj=1:n-1fori=j+1:nif(U(j,j)~=0)%判断主元元素是否不为0L(i,j)=U(i,j)/U(j,j);U(i,:)=U(i,:)-U(j,:)*U(i,j)/U(j,j);%U(j,j)为主元元素elsea=j+1;%令a等于j+1while((U(a,j)==0)&&(an))%判断主元元素所对的下一行元素是不是0，a是否小于na=a+1;%寻找下一个元素endtemp=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);endendendfork=1:nL(k,k)=1;%把L矩阵的对角线赋值为1endL%输出下三角矩阵LU%输出上三角矩阵UP%输出交换矩阵PA=inv(P)*L*Uelsedisp('error')endendif(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的第一列forj=2:size(A,2)m=zeros(size(A,1),1);%%取A的第一列fori=1:j-1R(i,j)=Q(:,i)'*A(:,j);%%q的转置乘以A的第j列向量m=m+R(i,j)*Q(:,i);%%q的转置乘以A的列向量endQ(:,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的列向量的单位化endQ%%输出正交矩阵QR%%输出上三角矩阵Rendif(x==3)%%判断如果x=3，那么将进行Householderreduction%%************Householderreduction***********%%disp('P*A=T')R=zeros(size(A,1));%%把R设置为矩阵维数等于矩阵的行数的全零方阵R1=zeros(size(A,1));%%把R1设为矩阵维数等于矩阵的行数的全零方阵M=A;%%将A赋给MP=eye(size(M,1));%%先将P矩阵设为维数等于M的单位矩阵fori=1:(size(M,1)-1)U=A;%%将A赋值给UU(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赋值给AA=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=[SV;FR];%%将SVFD合成矩阵R1elseR1=R;%%如果不满矩阵R的行数小于矩阵M的行数，则把R赋值给R1endP=R1*P;endP%%输出正交矩阵PT=P*M%%输出矩阵T,如果矩阵M的行数等于列数的话，T为上三角矩阵endif(x==4)%%判断x的值是否等于4，等于4则进行Givensreduction%%***********Givensreduction**********%%disp('P*A=R')U=A;%%将A赋值给Uw=size(A,1);%%w等于矩阵A的行数r=eye(w);%%将r设置为维数为w的单位矩阵fork=1:w-1m=eye(size(A,1));%%将m设置为维数等于A的行数单位矩阵fori=2:size(A,1)P=eye(size(A,1));a=0;%%将a是设置为0，方便求第一列前i个元素的平方和forj=1:iu=sqrt(a);a=a+A(j,1)^2;ends=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并赋值给mendA=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=zeros(size(m,1),w-size(m,1));p=[c,d;v,m];%%进行和并矩阵elsep=m;%%如果不满足矩阵m的行数小于w，则把m赋值给pendr=p*r;endP=r%%将r赋值给正交矩阵P，并输出PR=P*U%%输出矩阵R，若R的行数等于列数的话，R为上三角矩阵endend