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

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

实现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的行数