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

中科院矩阵分析与应用大作业实现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-1for 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%% 输出正交矩阵QR%% 输出上三角矩阵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赋给MP=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*