1、matlab上机作业报告计算初等反射阵用Householder变换法对矩阵A作正交分解连续函数最佳平方逼近等一、给定向量x0，计算初等反射阵Hk。1程序功能： 给定向量x0，计算初等反射阵Hk。2基本原理： 若的分量不全为零，则由确定的镜面反射阵H使得；当时，由有算法：(1)输入x，若x为零向量，则报错(2)将x规范化，如果M0，则报错同时转出停机否则(3)计算，如果，则(4)(5)计算(6)(7)(8)按要求输出,结束3变量说明：x 输入的n维向量；n n维向量x的维数；M M是向量x的无穷范数，即x中绝对值最大的一项的绝对值；p Householder初等变换阵的系数；u Househol
2、der初等变换阵的向量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 给定向量x0,计算Householder初等变换阵的p,u%程序功能：函数holder2给定向量x0,计算Householder初等变换阵的p,u；%输入：n维向量x；%输出：p,u。p是Householder初等变换阵的系数，% u是Householder初等变换阵的向量U。n=len
3、gth(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)n %如果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,
4、k:n)=eye(n-k+1)-pu*u; % 计算H(k:n,k:n)=I-pu*u；5使用示例:情形1:X为零向量 x=0,0,0,0; H=holderk(x,1)Error: M=0H = 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 = 检验: det(H)ans = H*xans = 情形4：(1)K值为3： x=4,3,2,1; H=holderk(x,3)H = 0 0 0 0 0 0 0 0 0 0
5、检验: det(H)ans = -1 H*xans = 0(2)K值为2: x=4,3,2,1; H=holderk(x,2)H = 0 0 0 0 0 0 det(H)ans = H*xans = 0二、设A为n阶矩阵，编写用Householder变换法对矩阵A作正交分解的程序。1程序功能：给定n阶矩阵A，通过本程序用Householder变换法对矩阵A作正交分解，得出AQR2基本原理： 任一实列满秩的mn矩阵A，可以分解成两个矩阵的乘积，即AQR，其中Q是具有法正交列向量的mn矩阵，R是非奇异的n阶上三角阵。算法：(1)输入n阶矩阵A(2)对，求Househoulder初等反射阵的。(3)
6、计算上三角阵R，仍然存储在A(4)计算正交阵Q (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阶矩阵
7、，% 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），对应初等反射阵的,u for jj=k:n %计算上三角阵R（仍存贮于A） t1=dot(u,A(k:n,jj)/p; %利用向量内积求和 A(k:n,jj)=A(k:n,jj)-t1*u; end for ii=1:n %计算正交矩阵Q t2=dot(u,Q(ii,k:n)/p
8、; %利用向量内积求和 Q(ii,k:n)=Q(ii,k:n)-t2*u; endend5使用示例:(1)A为3阶矩阵： A=1 2 3; 2 3 0; 3 4 5A = 1 2 3 2 3 0 3 4 5 q,r=qrhh(A)q = r = 0 检验： q*rans = (2)A为4阶矩阵： A=1 2 3 4; 2 3 0 1; 3 4 5 6;1 6 8 0A = 1 2 3 4 2 3 0 1 3 4 5 6 1 6 8 0 q,r=qrhh(A)q = r = 0 0 0 0 检验： q*rans = 0数值求解正方形域上的Poisson方程边值问题用MATLAB语言编写求解此辺值
9、问题的算法程序，采用下列三种方法，并比较三种方法的计算速度。1、用SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子；2、用块 Gauss-Sediel迭代法求解线性方程组Au=f；3、（预条件）共轭斜量法。用差分代替微分，对Poisson方程进行离散化，得到五点格式的线性方程组写成矩阵形式Au=f。其中 一 用SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子。1. 基本原理：Gauss-Seidel迭代法计算简单，但是在实际计算中，其迭代矩阵的谱半径常接近1，因此收敛很慢。为了克服这个缺点，引进一个加速因子（又称松弛因子）对Gauss-Seidel方法进行修正加速。假设
10、已经计算出第k步迭代的解（i=1，2，n），要求下一步迭代的解（i=1，2，n）。首先，用Gauss-Seidel迭代格式计算然后引入松弛因子，用松弛因子对和作一个线性组合。，i=1,2,n将二者合并成为一个统一的计算公式：2. 算法(1)Gauss-Seidel迭代法引入松弛因子w： 五点格式即为： (2)计算步骤： 第一步：给松弛因子赋初值w=,给场值u和场源b赋初值 第二步：用不同的w进行迭代计算。置error=0;计算在计算机上采用动态计算形式如果|du|error则error=|du|,如果errore,则停机，输出结果u,k.第三步：比较不同的w的迭代次数，用kk存放最小迭代次数，
11、用ww和uu存放相应的w及u。3 程序 u,k=SOR(u,b,w) %（被下面程序调用）%输入场初值u0、场源b及松弛因子w，通过五点差分格式进行迭代运算，%如果第k+1次的迭代结果与第k次的差小于精度，则可以近似认为第k+1次的迭代%结果是精确解，然后返回迭代次数k和迭代解function u,k=SOR(u,b,w) %输出迭代结果u，及迭代次数km=length(u); %m为u的维数h=1/(m-1); %h为步长N=10000;e=; %e为精度for k=1:N %k为记录迭代次数 error=0; for j=2:m-1 for jj=2:m-1 sum=4*u(jj,j)-u
12、(jj-1,j)-u(jj+1,j)-u(jj,j-1)-u(jj,j+1); du=w*(h2*b(jj,j)-sum)/4; %计算u的修正量 u(jj,j)= u(jj,j)+du; %修正u if errorabs(du), error=abs(du); end; %统计误差 end end if errork, kk=k;ww=w(i);uu=u;end %把最少迭代次数付给kk,及其w,u赋给ww,uuendt=toc %统计程序运算时间 4.计算结果： format short n=10; kk,ww,uu=SOR_5dianchafen(n) t = kk = 48ww = u
13、u = Columns 1 through 8 Columns 9 through 11 contourf (uu, DisplayName, uu); figure(gcf)图一 超松弛二 用块Jacobi迭代法求解线性方程组Au=f。1. 基本原理：对A做自然分解A=D-L-U=D-(L+U)其中D是有A的对角线元素组成的矩阵，L是由A的对角线以下元素组成的矩阵，U是由A得对角线以上元素组成的矩阵。于是将M=D,N=L+U,代入得到Dx=(L+U)x+b任取x的初值进行迭代2. 算法：（1）Gauss-Sediel迭代法原理五点差分格式： 因为A可以写成块状，即： 如果把每一条线上的节点看
14、作一个组，可以把Au=f表示成块状求解：（2）计算步骤： 第一步：给场值u和场源b赋初值，及定义 第二步：用公式，进行迭代计算 第三步：把第k次的u赋给ub，即ub=u；然后把第k+1次的u和ub进行比较，看是否达到精度，如果达到精度，则输出迭代次数k和精确解u。3. 程序k,u=kuai_GaussSeidel(n)%用块Gauss-Sediel迭代法求解正方形域上的Poisson方程边值问题%输入n,对x、y轴进行n等分；先确定场u的边界及场源b；%用k和u分别存放迭代次数和精确解function k,u=kuai_GaussSeidel(n) %对x、y轴进行n等分h=1/n; %步长u
15、=zeros(n+1,n+1); %对u赋初值u(1,1:n+1)=1;u(n+1,1:n+1)=1;u(1:n+1,1)=1;u(1:n+1,n+1)=1;b=zeros(n+1,n+1); %计算场源bfor i=2:n+1 for j=2:n+1 b(i,j)=(i-1)*(j-1)*h2; endend b=h2*b;for i=2:n b(2,i)=b(2,i)+u(1,i);b(i,n)=b(i,n)+u(i,n+1);b(n,i)= b(n,i)+u(n+1,i);b(i,2)= b(i,2)+u(i,1);endA=zeros(n-1,n-1); %定义矩阵的子块Afor i=
16、1:n-1 if i1, A(i,i-1)=-1; end if i2&jn, u(2:n,j)=pinv(A)*(u(2:n,j-1)+u(2:n,j+1)+b(2:n,j);end end error=max(max(abs(u-ub); %error是前后两次迭代结果的对应元素的最大误差 if error format short n=10; k,u=kuai_GaussSeidel(n)t = k = 93u = Columns 1 through 8 Columns 9 through 11 contourf (u, DisplayName, u); figure(gcf)图二 块G
17、auss-Sediel迭代法三 （预条件）共轭斜量法求解线性方程组Au=f。1.基本原理（1）预条件共轭斜量法原理（2）预优矩阵的选取 2. 算法3.程序%用J-PCG求解正方形域上的Poisson方程边值问题%输入n,对x、y轴进行n等分；先确定场u的边界及场源b；%用k和u分别返回迭代次数和精确解function k,u=J_CG(n)h=1/n; %h步长u=zeros(n+1,n+1); %对u赋初值u(1,1:n+1)=1;u(n+1,1:n+1)=1;u(1:n+1,1)=1;u(1:n+1,n+1)=1;b=zeros(n+1,n+1); %计算场源bfor i=2:n+1 fo
18、r j=2:n+1 b(i,j)=(i-1)*(j-1)*h2; endend b=h2*b;for i=2:n b(2,i)=b(2,i)+u(1,i);b(i,n)=b(i,n)+u(i,n+1); b(n,i)= b(n,i)+u(n+1,i);b(i,2)= b(i,2)+u(i,1);endz=zeros(n-1,n-1);r=zeros(n-1,n+1);p=zeros(n-1,n-1);bb=zeros(n-1,n-1);A=zeros(n-1,n-1); %定义矩阵的子块Afor i=1:n-1 if i1, A(i,i-1)=-1; end if i2&j1&i1&in-1,
19、r(:,i)= r(:,i)-a*(A*p(:,i)-p(:,i+1)-p(:,i-1);end if i=n-1,r(:,i)= r(:,i)-a*(A*p(:,i)-p(:,i-1);end end z(:,1:n-1)=pinv(A)*r(:,1:n-1); % kk=0;mm=0; for i=1:n-1 kk=kk+dot(z(:,i),r(:,i); mm=mm+dot(zz(:,i),rr(:,i); end beita=kk/mm; %beita是 p=z+beita*p; %对p进行更新，确定搜索方向 error=sqrt(norm(u-ub); %error是统计误差 if
20、 error format short n=10; k,u=J_CG(n)t = k = 37u = contourf (u, DisplayName, u); figure(gcf) 图三 J-PCG四 三种迭代法效率分析 有场的等值图可以看出，三种迭代方法的结果（达到相同的精度）比较一致，但是J-PCG只需要37次（耗时秒）迭代即可达到比较好的结果；超松弛迭代法需要48次（耗时秒）迭代即可达到满意的结果；块Gauss-Sediel迭代法需要93次（耗时）迭代才到达满意的结果。所以对此问题来说，超松弛迭代法和J-PCG法的效率要好于块Gauss-Sediel迭代法。一. 任务： 用MATLA
21、B语言编写连续函数最佳平方逼近的算法程序（函数式M文件）。并用此程序进行数值试验，写出计算实习报告。二. 程序功能要求： 在书中Page355或Page345的程序（见附一）的基础上进行修改，使其更加完善。要求算法程序可以适应不同的具体函数，具有一定的通用性。所编程序具有以下功能：1. 用Lengendre多项式做基，并适合于构造任意次数的最佳平方逼近多项式。可利用递推关系 2. 被逼近函数f(x)不用内联函数构造，而改用M文件建立数学函数。这样，此程序可通过修改建立数学函数的M文件以适用不同的被逼近函数（要学会用函数句柄）。3. 要考虑一般的情况。因此，程序中要有变量代换的功能。4. 计算组合系数时，计算函数的积分采用变步长复化梯形求积法（见附三）。5. 程序中应包括帮助文本和必要的注释语句。另外，程序中也要有必要的反馈信息。6. 程序输入：（
