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

中科院矩阵分析与应用大作业
实现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<n)) % 判断主元元素所对的下一行元素是不是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);
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)<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 %% 输出正交矩阵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)<w) %% 如果矩阵m的行数小于w
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。