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

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

实现LU分解 QR分解 Householder reduction、Givens reduction Matlab 代码：

function [] =juzhendazuoye A=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赋值给矩阵U L=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)&&(a

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); end end end for k=1:n

L(k,k)=1; % 把L矩阵的对角线赋值为1 end

L % 输出下三角矩阵L U % 输出上三角矩阵U P % 输出交换矩阵P A=inv(P)*L*U else 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); %% 把第一列赋值给a R(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 %% 输出上三角矩阵R end

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)

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 %% 输出正交矩阵P

T=P*M %% 输出矩阵T,如果矩阵M的行数等于列数的话，T为上三角矩阵 end

if(x==4) %% 判断x的值是否等于4，等于4则进行Givens reduction

%%***********Givens reduction**********%% disp('P*A=R')

U=A; %% 将A赋值给U

w=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 end

A=m*A; %% 矩阵m*A赋值给A A=A(2:size(A,1),2:size(A,2)); %% 取A的子矩阵 if(size(m,1)

c=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]; %% 进行和并矩阵 else

p=m; %% 如果不满足矩阵m的行数小于w，则把m赋值给p end

r=p*r; end

P=r %% 将r赋值给正交矩阵P，并输出P

R=P*U %% 输出矩阵R，若R的行数等于列数的话，R为上三角矩阵 end end

百度搜索“77cn”或“免费范文网”即可找到本站免费阅读全部范文。收藏本站方便下次阅读，免费范文网，提供经典小说综合文库中科院矩阵分析与应用大作业在线全文阅读。