一、给定向量x≠0，计算初等反射阵Hk。 1．程序功能：
给定向量x≠0，计算初等反射阵Hk。 2．基本原理：
若x??x,x,?,x??R的分量不全为零，则由
???sign(x1)x2??u?x??e1?(x1??,x2,?,xn)?1?2 ??u??(??x1)?22?T?uu?1T?I??uu?H?I?22u2??确定的镜面反射阵H使得Hx???e?y；当(1?k?n)时，由
n?21/2??sign(x)((x))?kki?i?k??u(k)?(0,?,0,xk??k,xk?1,?,xn)T?Rn ????1(u(k))Tu(k)??(??x)=?(u(k))Tkkkk?k2??1(k)(k)T?Hk?I??ku(u)Tn有Hkx?(x1,x2,?,xk?1,??k,0,?,0)?R
算法：
(1)输入x，若x为零向量，则报错 (2)将x规范化，
M?max?x,x,?,x?
如果M＝0，则报错同时转出停机 否则xi?xiM,i?1,2,?,n (3)计算??x2，如果x1?0，则????
(4)???(??x1) (5)计算u?x,u(1)?x1?? (6)H?I??uu (7)y?(?M?,0,?,0)
?1T
(8)按要求输出,结束
3．变量说明：
x － 输入的n维向量；
n － n维向量x的维数；
M － M是向量x的无穷范数，即x中绝对值最大的一项的绝对值； p － Householder初等变换阵的系数ρ； u － Householder初等变换阵的向量U s x n
－ 向量x的二范数； － 输入的n维向量； － n维向量x的维数；
p － Householder初等变换阵的系数ρ； u － Householder初等变换阵的向量U
k － 数k，H*x=y，使得y的第k+1项到最后项全为零； 4．程序代码：
(1)
function [p,u]=holder2(x)
%HOLDER2 给定向量x≠0,计算Householder初等变换阵的p,u
%程序功能：函数holder2给定向量x≠0,计算Householder初等变换阵的p,u； %输入：n维向量x；
%输出：[p,u]。p是Householder初等变换阵的系数ρ， % u是Householder初等变换阵的向量U。 n=length(x); % 得到n维向量x的维数； p=1;u=0; % 初始化p,u； M=max(abs(x)); if M==0
% 得到向量x的无穷范数，即x中绝对值最大的一项的绝对值；
% 如果x=0，提示出错,程序终止;
disp('Error: M=0');
return; else
x=x/M; % 规范化 end;
s=norm(x); if x(1)<0 s=-s; end u=x;
% 求x的二范数
% 首项为负，s值要变号
% 除首项外，其余各项x,u相同 % 计算u的首项 % 计算p
u(1)=s+x(1); p=s*u(1); (2)
if n==1 u=0; end % 若x是1×1维向量，则u=0
function H=holderk(x,k)
%HOLDERK 给定向量x≠0,数k，计算初等反射阵Hk，使HkX=Y，其中Y的第k+1项到
最后项全为零；
%程序功能：函数holderk给定向量x≠0，数k，计算初等反射阵Hk，使HkX=Y，%程序功能：函数holder2给定向量x≠0,计算Householder初等变换阵的p,u； %输入：n维向量x,数k；
%输出：H。H是Householder初等变换阵，H*x=y，使得y的第k+1项到最后项全为零； %引用函数：holder2；
n=length(x); % 得到n维向量x的维数； if k>n %如果k值溢出，报错； disp('Error: k>n');
end
H=eye(n); % 初始化H，并使H(1:k,1:k)=I；
[p,u]=holder2(x(k:n)); % 得到计算Householde初等变换阵的系数ρ、向量U； H(k:n,k:n)=eye(n-k+1)-p\%u*u'; % 计算H(k:n,k:n)=I-p\%u*u'； 5．使用示例: 情形1: X为零向量 >> x=[0,0,0,0]'; >> H=holderk(x,1) Error: M=0 H =
1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 情形2： K值溢出： >> x=[1,2,3,4]';
>> H=holderk(x,5) Error: k>n 情形3： K值为1： >> x=[2,3,4,5]'; >> H=holderk(x,1) H =
-0.2722 -0.4082 -0.5443 -0.6804 -0.4082 0.8690 -0.1747 -0.2184 -0.5443 -0.1747 0.7671 -0.2911 -0.6804 -0.2184 -0.2911 0.6361 检验: >> det(H)
ans =
-1.0000 >> H*x ans =
-7.3485 0.0000 0.0000 0.0000 情形4： (1)K值为3： >> x=[4,3,2,1]'; >> H=holderk(x,3) H =
1.0000 0 0 1.0000 0 0 0 0 检验:
>> det(H)
ans = -1
>> H*x
ans =
4.0000 3.0000 -2.2361 0 (2)K值为2: >> x=[4,3,2,1]'; >> H=holderk(x,2)
0 0 0 0 -0.8944 -0.4472 -0.4472 0.8944
H =
1.0000 0 0 0 0 -0.8018 -0.5345 -0.2673 0 -0.5345 0.8414 -0.0793 0 -0.2673 -0.0793 0.9604
>> det(H) ans =
-1.0000 >> H*x
ans =
4.0000 -3.7417 0.0000
0
二、设A为n阶矩阵，编写用Householder变换法对矩阵A作正交分解的程序。 1．程序功能：
给定n阶矩阵A，通过本程序用Householder变换法对矩阵A作正交分解，得出A＝QR 2．基本原理：
任一实列满秩的m×n矩阵A，可以分解成两个矩阵的乘积，即A＝QR，其中Q是具有法正交列向量的m×n矩阵，R是非奇异的n阶上三角阵。 算法：
(1)输入n阶矩阵A
(2)对A?k:n,k?，求Househoulder初等反射阵的?k,u(3)计算上三角阵R，仍然存储在A
HkA(k)(k)。k?1,2,?,n?1
?A(k?1)
j?k,k?1,?,n
1?n?k?(k)?tj???uiaij??k?i?k? aij(k?1)?aij(k)?tjui?k?(4)计算正交阵Q