数值分析matlab上机作业报告.docVIP

一、给定向量x≠0，计算初等反射阵Hk。 1．程序功能： 给定向量x≠0，计算初等反射阵Hk。 2．基本原理： 若的分量不全为零，则由 确定的镜面反射阵H使得；当时，由 有 算法： (1)输入x，若x为零向量，则报错 (2)将x规范化， 如果M＝0，则报错同时转出停机 否则 (3)计算，如果，则 (4) (5)计算 (6) (7) (8)按要求输出,结束 3．变量说明： x － 输入的n维向量； n － n维向量x的维数； M － M是向量x的无穷范数，即x中绝对值最大的一项的绝对值； p － Householder初等变换阵的系数ρ； u － Householder初等变换阵的向量U s － 向量x的二范数； x － 输入的n维向量； 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)); % 得到向量x的无穷范数，即x中绝对值最大的一项的绝对值； if M==0 % 如果x=0，提示出错,程序终止; disp(Error: M=0); return; else x=x/M; % 规范化 end; s=norm(x); % 求x的二范数 if x(1)0 % 首项为负，s值要变号 s=-s; end u=x; % 除首项外，其余各项x,u相同 u(1)=s+x(1); % 计算u的首项 p=s*u(1); % 计算p if n==1 u=0; end % 若x是1×1维向量，则u=0 (2) 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 kn %如果k值溢出，报错； disp(Error: kn); 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: kn 情形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