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

一、给定向量x≠0，计算初等反射阵Hk。1．程序功能：给定向量x≠0，计算初等反射阵Hk。2．基本原理：若xxxRx,,,的分量不全为零，则由12112212122()x(,,,)1()22nTTsignxexxxxuxuuuHIIuuu确定的镜面反射阵H使得yeHx；当(1)kn时，由21/2k()T1()()()k1()()()(())(0,,0,,,,)1()()=()2()nkiikknkkknkTkkTkkkkkkTkksignxxxxxxuRuuuHIuu有T121(,,,,,0,,0)nkkkxxxHxR算法：(1)输入x，若x为零向量，则报错(2)将x规范化，xxxM,,,max如果M＝0，则报错同时转出停机否则niMxxii,,2,1,(3)计算2x，如果01x，则(4))(1x(5)计算1,(1)xuxu(6)1THIuu(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中绝对值最大的一项的绝对值；ifM==0%如果x=0，提示出错,程序终止;disp('Error:M=0');return;elsex=x/M;%规范化end;s=norm(x);%求x的二范数ifx(1)0%首项为负，s值要变号s=-s;endu=x;%除首项外，其余各项x,u相同u(1)=s+x(1);%计算u的首项p=s*u(1);%计算pifn==1u=0;end%若x是1×1维向量，则u=0(2)functionH=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的维数；ifkn%如果k值溢出，报错；disp('Error:kn');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=1000010000100001情形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.40820.8690-0.1747-0.2184-0.5443-0.17470.7671-0.2911-0.6804-0.2184-0.29110.6361检验:det(H)ans=-1.0000H*xans=-7.34850.00000.00000.0000情形4：(1)K值为3：x=[4,3,2,1]';H=holderk(x,3)H=1.000000001.00000000-0.8944-0.447200-0.44720.8944检验:det(H)ans=-1H*xans=4.00003.0000-2.23610(2)K值为2:x=[4,3,2,1]';H=holderk(x,2)H=1.00000000-0.8018-0.5345-0.26730-0.53450.8414-0.07930-0.2673-0.07930.9604det(H)ans=-1.0000H*xans=4.0000-3.74170.00000二、设A为n阶矩阵，编写用Householder变换法对矩阵A作正交分解的程序。1．程序功能：给定n阶矩阵A，通过本程序用Householder变换法对矩阵A作正交分解，得出A＝QR2．基本原理：任一实列满秩的m×n矩阵A，可以分解成两个矩阵的乘积，即A＝QR，其中Q是具有法正交列向量的m×n矩阵，R是非奇异的n阶上三角阵。算法：(1)输入n阶矩阵A(2)对:,Aknk，求Househoulder初等反射阵的)(,kku。1,,2,1nk(3)计算上三角阵R，仍然存储在A)1()(kkkAAHnkkj,,1,kijkijkijnkikijkikjutaaaut)()1()(1(4)计算正交阵Q)1()()1(kkkQHQIQni,,2,1kjikijkijnkjkjkijkiutqqnkkjuqt)()1()(,,1,1(5)按要求输出,结束3．变量说明：A－输入的n阶矩阵，同时用于存储上三角阵R；n－矩（方）阵A的阶数；Q－Q是具有法正交列向量的n阶矩阵；p,u－向量A（k:n,k），对应初等反射阵的ρ,uk，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；fork=1:n-1[p,u]=holder2(A(k:n,k));%向量A（k:n,k），对应初等反射阵的ρ,uforjj=k:n%计算上三角阵R（仍存贮于A）t1=dot(u,A(k:n,jj))/p;%利用向量内积求和A(k:n,jj)=A(k:n,jj)-t1*u;endforii=1:n%计算正交矩阵Qt2=dot(u,Q(ii,k:n))/p;%利用向量内积求和Q(ii,k:n)=Q(ii,k:n)-t2*u';endend5．使用示例:(1)A为3阶矩阵：A=[123;230;345]A=123230345[q,r]=qrhh(A)q=-0.26730.87290.4082-0.53450.2182-0.8165-0.8018-0.43640.4082r=-3.7417-5.3452-4.810700.65470.4364-0.00000.00003.2660检验：q*rans=1.00002.00003.00002.00003.00000.00003.00004.00005.0000(2)A为4阶矩阵：A=[1234;2301;3456;1680]A=1234230134561680[q,r]=qrhh(A)q=-0.25820.0597-0.2660-0.9268-0.5164-0.10450.8434-0.1049-0.7746-0.2688-0.46620.3323-0.25820.9556-0.02220.1399r=-3.8730-6.7132-6.7132-6.196804.46476.4805-1.47830-0.0000-3.3070-3.017800.00000-1.8187检验：q*rans=1.00002.00003.00004.00002.00003.0000-0.00001.00003.00004.00005.00006.00001.00006.00008.00000数值求解正方形域上的Poisson方程边值问题用MATLAB语言编写求解此辺值问题的算法程序，采用下列三种方法，并比较三种方法的计算速度。1、用SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子；2、用块Gauss-Sediel迭代法求解线性方程组Au=f；3、（预条件）共轭斜量法。用差分代替微分，对Poisson方程进行离散化，得到五点格式的线性方程组2,1,1,,1,11,1,,1,142,11,1ijijijijijijjNjiiNuuuuuhfijNuuuuijN，写成矩阵形式Au=f。其中2222(,),0,1(0,)(1,)(,0)(,1)1uufxyxyxyxyuyuyuxux2233NNvbvbufvb4114114iiA22,23,2,232,33,3,32,3,,222,23,2,22,11,23,1,11,2232,33,3,31,31,3(,,...,),(,,...,),......,(,,...,)(,,...,)(,,...,),(,,...,)(,0,...,),......,TTNNTNNNNNTTNNNTTNNNvuuuvuuuvuuubhfffuuuuubhfffufb22,3,,2,11,3,1,11,,(,,...,)(,,...,)1,10,0.1,,2,3,...,TTNNNNNNNNNNNijijhfffuuuuuhNhNfxyijN取则212331,1NNAIIAAIIA一用SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子。1.基本原理：Gauss-Seidel迭代法计算简单，但是在实际计算中，其迭代矩阵的谱半径常接近1，因此收敛很慢。为了克服这个缺点，引进一个加速因子（又称松弛因子）对Gauss-Seidel方法进行修正加速。假设已经计算出第k步迭代的解（i=1，2，···，n），要求下一步迭代的解（i=1，2，···，n）。首先，用Gauss-Seidel迭代格式计算)1(kx然后引入松弛因子，用松弛因子对和作一个线性组合。)()1()1()1(kikikixxx，i=1,2,…,n将二者合并成为一个统一的计算公式：2.算法(1)Gauss-Seidel迭代法引入松弛因子w：五点格式即为：(2)计算步骤：第一步：给松弛因子赋初值w=1.1～1.8,给场值u和场源b赋初