matlab上机作业报告(计算初等反射阵-用Householder变换法对矩阵A作正交分解-连续函数最佳平方逼近等)

一、给定向量x ≠0，计算初等反射阵H k 。

1．程序功能：给定向量x ≠0，计算初等反射阵H k 。

2．基本原理： 若()xx x R x ∈=,,, 的分量不全为零，则由12112212122()x (,,,)1()22n T T sign x e x x x x σσσρσσρ-⎧=⎪=+=+⎪⎪⎪==+⎨⎪⎪=-=-⎪⎪⎩u x u uu H I I uuu 确定的镜面反射阵H 使得y eHx =-=σ；当(1)k n ≤<时，由21/2k ()T 1()()()k 1()()()(())(0,,0,,,,)1()()=()2()nk i i kk nk k k n k T k k Tk k k kk k T k k sign x x x x x x σσρσσσρ=+-⎧=⎪⎪⎪=+∈⎨⎪==+⎪⎪=-⎩∑u R u u u H I u u 有T 121(,,,,,0,,0)n k k k x x x σ-=-∈H x R算法：(1)输入x ，若x 为零向量，则报错 (2)将x 规范化，{}x x x M ,,,m ax =如果M ＝0，则报错同时转出停机 否则n i M x x i i ,,2,1, =←(3)计算2x =σ，如果0<1x ，则σσ-= (4))(1x +=σσρ (5)计算1,(1)x σ==+u x u (6)1Tρ-=-H I uu (7)(M ,0,,0)σ=-y(8)按要求输出,结束3．变量说明：x －输入的n维向量；n －n维向量x的维数；M －M是向量x的无穷范数，即x中绝对值最大的一项的绝对值；p －Householder初等变换阵的系数ρ；u －Householder初等变换阵的向量Us －向量x的二范数；x －输入的n维向量；n －n维向量x的维数；p －Householder初等变换阵的系数ρ；u －Householder初等变换阵的向量Uk －数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;elsex=x/M; % 规范化end;s=norm(x); % 求x的二范数if x(1)<0 % 首项为负，s值要变号s=-s;endu=x; % 除首项外，其余各项x,u相同u(1)=s+x(1); % 计算u的首项p=s*u(1); % 计算pif 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 k>n %如果k值溢出，报错；disp('Error: k>n');endH=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=0H =1 0 0 00 1 0 00 0 1 00 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 =检验:>> det(H)>> H*xans =情形4：(1)K值为3：>> x=[4,3,2,1]';>> H=holderk(x,3)H =0 0 0 0 0 0 0 00 0检验:>> det(H)ans =-1>> H*xans =(2)K值为2:>> x=[4,3,2,1]';>> H=holderk(x,2)0 0 0 0 0 0>> det(H)ans =>> H*xans =二、设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 k u ρ。

1,,2,1-=n k(3)计算上三角阵R ，仍然存储在A)1()(+→k k k A A Hn k k j ,,1, +=()()k i j k ij k ij n k i k ij k i k j u t a a a u t -=⎪⎭⎫ ⎝⎛=+=∑)()1()(1ρ (4)计算正交阵Q)1()()1(+→=k k k QH Q I Qn i ,,2,1 =()()k j i k ij k ij nkj k j k ijki u t q q nk k j u qt -=+==+=∑)()1()(,,1,1ρ(5)按要求输出,结束3．变量说明：A － 输入的n 阶矩阵，同时用于存储上三角阵R ； n － 矩（方）阵A 的阶数；Q － Q 是具有法正交列向量的n 阶矩阵； p,u － 向量A （k:n,k ），对应初等反射阵的ρ,u k ，jj ，ii － 循环变量；t1 － 计算上三角阵R 的系数tj ； t2 － 计算正交矩阵Q 的系数ti ； 4．程序代码：function [Q,A]=qrhh(A)%QRHH 用Householder 变换法对n 阶矩阵A 作正交分解A=QR ；%程序功能：函数qrhh 用Householder 变换法对矩阵A 作正交分解A=QR ； %输入：n 阶矩阵A ；%输出：[Q,A]。

Q 是具有法正交列向量的n 阶矩阵，% A （即R ）是非奇异的n 阶上三角阵，仍用输入的矩阵A 存储。

%引用函数：% holder2；示例 [p,u]=holder2(x); [n,n]=size(A); %求矩（方）阵A 的阶数； Q=eye(n); %构造正交矩阵Q(1)=I ； for k=1:n-1[p,u]=holder2(A(k:n,k)); %向量A （k:n,k ），对应初等反射阵的ρ,ufor jj=k:n %计算上三角阵R （仍存贮于A ） t1=dot(u,A(k:n,jj))/p; %利用向量内积求和 A(k:n,jj)=A(k:n,jj)-t1*u; endfor ii=1:n %计算正交矩阵Qt2=dot(u,Q(ii,k:n))/p; %利用向量内积求和 Q(ii,k:n)=Q(ii,k:n)-t2*u'; end end5．使用示例: (1)A 为3阶矩阵：>> A=[1 2 3; 2 3 0; 3 4 5]A =1 2 32 3 03 4 5>> [q,r]=qrhh(A)q =r =检验：>> q*rans =(2)A为4阶矩阵：>> A=[1 2 3 4; 2 3 0 1; 3 4 5 6;1 6 8 0]A =1 2 3 42 3 0 13 4 5 61 6 8 0>> [q,r]=qrhh(A)q =r =0 00 0 检验： >> q*rans =数值求解正方形域上的Poisson 方程边值问题2222(,),0,1(0,)(1,)(,0)(,1)1u u f x y xy x y x y u y u y u x u x ⎧⎛⎫∂∂-+==<<⎪ ⎪∂∂⎨⎝⎭⎪====⎩用MATLAB 语言编写求解此辺值问题的算法程序，采用下列三种方法，并比较三种方法的计算速度。

1、用SOR 迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子；2、用块 Gauss-Sediel 迭代法求解线性方程组Au=f ；3、（预条件）共轭斜量法。

用差分代替微分，对Poisson 方程进行离散化，得到五点格式的线性方程组2,1,1,,1,11,1,,1,142,11,1i j i j i j i j i j ij j N ji i N u u u u u h f i j N u u u u i j N -+-+++----=≤≤====≤≤+，写成矩阵形式Au=f 。

其中一 用SOR 迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子。

1. 基本原理：Gauss-Seidel 迭代法计算简单，但是在实际计算中，其迭代矩阵的谱半径常接近1，因此收敛很慢。

为了克服这个缺点，引进一个加速因子（又称松弛因子）对Gauss-Seidel2233N N v b v b u f v b ⎛⎫⎛⎫⎪ ⎪ ⎪⎪== ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭4114114ii A -⎛⎫ ⎪- ⎪= ⎪- ⎪-⎝⎭22,23,2,232,33,3,32,3,,222,23,2,22,11,23,1,11,2232,33,3,31,31,3(,,...,),(,,...,),......,(,,...,)(,,...,)(,,...,),(,,...,)(,0,...,),......,T T N N TN N N N N T T N N N T T N N N v u u u v u u u v u u u b h f f f u u u u u b h f f f u f b ++====+++=+=22,3,,2,11,3,1,11,,(,,...,)(,,...,)1,10,0.1,,2,3,...,T TN N N N N N N N N N N i j i j h f f f u u u u u h N h N f x y i j N+++++++=====取则212331,1N N A I I A A I IA ---⎛⎫ ⎪-⎪= ⎪- ⎪⎪-⎝⎭方法进行修正加速。