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

中科院矩阵分析与应用大作业
实现LU分解QR分解Householder reduction、Givens reduction
Matlab 代码：
function [] =j uzhendazuoye
A=input ('请输入•个矩阵A=');
2 Gram-Schmidt 分解
3 Householder reduction 4
x=input (*请输入序号1 LU分解
Givens reduction: 1);
if (>:==!)
壮mmm分解mm%%
disp('PA=LU1)
m=size(A,1); %nt等于矩阵A的行数
n=size(A,2); %n等于矩阵A的列数
if (m==n) % 刊斯NA是不足方阵
% 如果矩阵A不是方阵那么就输出"error"
U=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 z j)=U(i z j) /U(j z j)；
U(i f :)=U(i, :)-U(j, j)/U(j z j); % U(j, j)为主元元素
else
a=j+l；% 令 a 等于j + 1
while ( (U (a, j ) ==0) && (a<n) ) %判断主元元素所刘•的卜行九索, 上0, a是否小于n
a=a+l; %寻找下•个元素
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(jr :)；
P(j,:)=P(a,:); %交换矩阵的两行交换
P(a,：)=q；
L(i z j)=U(i z j)/U(j r j)；
U(i, :)=U(i, :)-U(j, :)*U(i f j)/U(j r j);
end
end
end for k=l:n
L(k,k)=l; %把L矩阵的对角线赋值为1
end
L %输出下三角矩阵
u %输岀上三角矩阵u
p %输出交换矩阵p
A=inv(P)*L*U
else disp('error 1)
end
end
if (x==2) %%判断如果那么将执行sohmid分解%%**************Gram-Schmidc 正交分解mm
disp(,A=Q*R,)
Q=zeros (size (A, 1) , size (A z 2) ) ; %% 先把Q 设为全零矩阵
R=zeros(size(A,2)); %% R设置为全零矩阵
a=A(：,1)；%% 把第•列赋值给a
R(1, 1)=norm(a); %% 求第•列列向量的模值
a=a/norm(a); %% 求第•列列向量的单位向量
Q(：,l)=a; %% 把a賦值给Q的第•列
for j=2:size(A z2)
m=zeros (size (A, 1),1); %% 取 A 的第•列
for i=l:j-1
R(i, j)=Q(：zi) A(:, j); %% q的转置乘以A的第j列向量
m=m+R(i z j ) *Q ( : z i); %% q的转置乘以A的列向量
end
Q ( : , j ) =A (:，j ) -m; %%企的第：列减去q(i)和A (: f j )的内积和
R( jr j) =norm (Q ( ：z 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 f)
R=zeros(size(A,1)); %%把R设置为矩阵维数等于矩阵的行数的全零方阵
Rl = zeros (size(A,1)); %%把R1设为矩阵维数等于矩阵的行数的全零方阵
M=A; %%将A赋给M
P=eye(size(M,1)); %%先将E矩阵设为维数等于M的单位矩阵for 1=1:(size(M z1)-1)
U=A; %%将A賦值给U
U(l z 1)=U(1Z 1) -norm (U(:, 1)); 锂将U的第•列的第彳了元素减去U的; ］列
向量的模值
R=eye(size(U,1))-2*U(:,1)*U(:,1)f/(U(:z l) ** U(:z l))；%%
I-2*U(:,1)*U(:,1) f/(U(:, 1) U(:z 1)
A=R*A;
%% R乘以入賦值给A
A=A (2 : size (A, 1) , 2 : size (A, 2) ) ; %% 取g 的 /矩阵
if (size (R, 1) Size (M, 1) ) %%判断矩阵R的行数是否小于矩阵M的行数.
如果小于进行下步：
S=eye (size (M, l)-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) z size (M, 1) -size (R z 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的 f J•列数的话，T为上
三角矩阵
end
if (x==4) “判断瓦的值是笛需于4,等于4则进行Givzs
reduction
***★*★★***Givens reduction**********%%
disp(,P*A=R,)
w=size(A,1); %% w等于矩阵A的行数
U=A; %% 将A赋值给U
r=eye(w)；%% 将r设置为维数为w的单位矩阵
for k=l:w-1
m=eye(size(A,1)); %% 将m设置为维数等于A的行数单位矩阵for i=2:size(A z l)
P=eye(size(A,1));
a=0; $$将&是设置为0,方便求第•列前i个元素的平方和for j=l:i u=sqrt(a); a=a+A(j z l)A2;
end
s=sqrt(a)；%%将第-列前i个元素的平方开根
P(l z l)=u/s; %%将u/s賦值给旋转矩阵P的第•行的第•列
P(i,i)=u/s; %%将U/S賦值给旋转矩阵P的第i行和第i列
P(i z l)=-A(i z l)/s；%%将-A(i,D赋值非P的第i行的第•列
P(l z i)=A(i,l)/s；%%将A(i z i)赋值给P的第•行的第i列
m=P*m; %% P乘以矩阵m并賦值给m
end
A=m*A; %%矩阵赋值给A
A=A (2 : size (A, 1),2: size (A, 2) ) ; %% A 的/妙阵
if (size(m z1)<w) %%如果矩阵in的行数小于w
c=eye(w-size(m,1)); 将o设置为维数等于w-矩阵m的行数的单位矩阵d=zeros(w-size(m,1),size(m,1));
v=zeros(size(m,1),w-size(m,1));
w=size(A,1); %% w等于矩阵A的行数
P= [c, d; v z m]; %%进行和并矩阵else
p=r end
r=p*r; end
P=r
R=P*U end %%如果不满足矩阵m的行数小于w,则把%%将r赋值给正交矩阵P,并输出P
end
赋值给p %%输出矩阵R，若R的行•数等于列数的话，R为上三角矩阵。