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值溢出，报错； knH=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）-pu*u; % 计算H（k:n）=I-pu*u5使用
2、示例:情形1:X为零向量 x=0,0,0,0 H=holderk（x,1） 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）n情形3：K值为1： x=2,3,4,5 -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 -7.3485 0.0000情形4
3、：（1）K值为3： x=4,3,2,1 H=holderk（x,3） 1.0000 0 0 0 0 1.0000 0 0 0 0 -0.8944 -0.4472 0 0 -0.4472 0.8944 -1 4.0000 3.0000 -2.2361 0（2）K值为2: H=holderk（x,2） 0 -0.8018 -0.5345 -0.2673 0 -0.5345 0.8414 -0.0793 0 -0.2673 -0.0793 0.9604 -3.7417二、设A为n阶矩阵，编写用Householder变换法对矩阵A作正交分解的程序。给定n阶矩阵A，通过本程序用Householder变换
4、法对矩阵A作正交分解，得出AQR 任一实列满秩的mn矩阵A，可以分解成两个矩阵的乘积，即AQR，其中Q是具有法正交列向量的mn矩阵，R是非奇异的n阶上三角阵。（1）输入n阶矩阵A（2）对，求Househoulder初等反射阵的。（3）计算上三角阵R，仍然存储在A（4）计算正交阵Q （5）按要求输出,结束A 输入的n阶矩阵，同时用于存储上三角阵R；n 矩（方）阵A的阶数；Q Q是具有法正交列向量的n阶矩阵；p,u 向量A（k:n,k），对应初等反射阵的,uk，jj，ii 循环变量；t1 计算上三角阵R的系数tj；t2 计算正交矩阵Q的系数ti；function Q,A=qrhh（A）%QRHH
5、用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: for jj=k:n %计算上三角阵R（仍存贮于A） t1=dot（u,A（k:n,jj）/p; %利用向量内积求和
6、A（k:n,jj）=A（k:n,jj）-t1*u; end for ii=1:n %计算正交矩阵Q t2=dot（u,Q（ii,k:n）/p; %利用向量内积求和 Q（ii,k:n）=Q（ii,k:n）-t2*u（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 = -0.2673 0.8729 0.4082 -0.5345 0.2182 -0.8165 -0.8018 -0.4364 0.4082r = -3.7417 -5.3452 -4.8107 0 0.6547 0.4364 -0.0000 0.0000
7、 3.2660检验： q*r 1.0000 2.0000 3.0000 2.0000 3.0000 0.00003.0000 4.0000 5.0000（2）A为4阶矩阵： A=1 2 3 4; 2 3 0 1; 3 4 5 6;1 6 8 0 1 2 3 4 2 3 0 1 3 4 5 6 1 6 8 0 -0.2582 0.0597 -0.2660 -0.9268 -0.5164 -0.1045 0.8434 -0.1049 -0.7746 -0.2688 -0.4662 0.3323 -0.2582 0.9556 -0.0222 0.1399 -3.8730 -6.7132 -6.713
8、2 -6.1968 0 4.4647 6.4805 -1.4783 0 -0.0000 -3.3070 -3.0178 0 0.0000 0 -1.8187 1.0000 2.0000 3.0000 4.0000 2.0000 3.0000 -0.0000 1.0000 3.0000 4.0000 5.0000 6.0000 1.0000 6.0000 8.0000 0数值求解正方形域上的Poisson方程边值问题用MATLAB语言编写求解此辺值问题的算法程序，采用下列三种方法，并比较三种方法的计算速度。1、用SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子；2、用块 Gauss-
9、Sediel迭代法求解线性方程组Au=f；3、（预条件）共轭斜量法。用差分代替微分，对Poisson方程进行离散化，得到五点格式的线性方程组写成矩阵形式Au=f。其中一 用SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子。1. 基本原理：Gauss-Seidel迭代法计算简单，但是在实际计算中，其迭代矩阵的谱半径常接近1，因此收敛很慢。为了克服这个缺点，引进一个加速因子（又称松弛因子）对Gauss-Seidel方法进行修正加速。假设已经计算出第k步迭代的解（i=1，2，n），要求下一步迭代的解（i=1，2，n）。首先，用Gauss-Seidel迭代格式计算然后引入松弛因子，用松弛因
10、子对和作一个线性组合。，i=1,2,n将二者合并成为一个统一的计算公式：2.算法（1）Gauss-Seidel迭代法引入松弛因子w：五点格式即为：（2）计算步骤： 第一步：给松弛因子赋初值w=1.11.8,给场值u和场源b赋初值 第二步：用不同的w进行迭代计算。置error=0;计算在计算机上采用动态计算形式如果|du|error则error=|du|,如果errore,则停机，输出结果u,k.第三步：比较不同的w的迭代次数，用kk存放最小迭代次数，用ww和uu存放相应的w及u。3 程序 u,k=SOR（u,b,w） %（被下面程序调用）%输入场初值u0、场源b及松弛因子w，通过五点差分格式进
11、行迭代运算，%如果第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=0.0000001; %e为精度N %k为记录迭代次数 error=0; for j=2:m-1 for jj=2:m-1 sum=4*u（jj,j）-u（jj-1,j）-u（jj+1,j）-u（jj,j-1）-u（jj,j+1）; du=w*（h2*b（jj,j）-sum）/4; %计算u的修正
12、量 u（jj,j）= u（jj,j）+du; %修正u if errork, kk=k;ww=w（i）;uu=u;end %把最少迭代次数付给kk,及其w,u赋给ww,uut=toc %统计程序运算时间4.计算结果： format short n=10; kk,ww,uu=SOR_5dianchafen（n） t = 0.0310kk = 48ww = 1.6000uu = Columns 1 through 8 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0011 1.0022 1.0031 1.0039
13、1.0044 1.0047 1.0045 1.0000 1.0022 1.0042 1.0061 1.0076 1.0087 1.0091 1.0088 1.0000 1.0031 1.0061 1.0088 1.0110 1.0126 1.0133 1.0128 1.0000 1.0039 1.0076 1.0110 1.0138 1.0159 1.0168 1.0162 1.0000 1.0044 1.0087 1.0126 1.0159 1.0183 1.0194 1.0189 1.0000 1.0047 1.0091 1.0133 1.0168 1.0194 1.0208 1.0203
14、 1.0000 1.0045 1.0088 1.0128 1.0162 1.0189 1.0203 1.0201 1.0000 1.0037 1.0073 1.0107 1.0136 1.0160 1.0174 1.0175 1.0000 1.0023 1.0045 1.0066 1.0084 1.0100 1.0110 1.0113 Columns 9 through 11 1.0000 1.0000 1.0000 1.0037 1.0023 1.0000 1.0073 1.0045 1.0000 1.0107 1.0066 1.0000 1.0136 1.0084 1.0000 1.016
15、0 1.0100 1.0000 1.0174 1.0110 1.0000 1.0175 1.0113 1.0000 1.0155 1.0103 1.0000 1.0103 1.0072 1.0000 contourf （uu, DisplayName, uu figure（gcf）图一 超松弛二 用块Jacobi迭代法求解线性方程组Au=f。对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
16、-Sediel迭代法原理五点差分格式：因为A可以写成块状，即： 如果把每一条线上的节点看作一个组，可以把Au=f表示成块状求解：给场值u和场源b赋初值，及定义用公式，进行迭代计算 第三步：把第k次的u赋给ub，即ub=u；然后把第k+1次的u和ub进行比较，看是否达到精度，如果达到精度，则输出迭代次数k和精确解u。3.程序k,u=kuai_GaussSeidel（n）%用块Gauss-Sediel迭代法求解正方形域上的Poisson方程边值问题先确定场u的边界及场源b；%用k和u分别存放迭代次数和精确解function k,u=kuai_GaussSeidel（n） %对x、y轴进行n等分 %
17、步长u=zeros（n+1,n+1）;b=h2*b; 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）;A=zeros（n-1,n-1）; %定义矩阵的子块A if i1, A（i,i-1）=-1; if i2&jn, u（2:n,j）=pinv（A）*（u（2:n,j-1）+u（2:n,j+1）+b（2:n,j）; error=max（max（abs（u-ub）; %error是前后两次迭代结果的对应元素的最大误差=e, break; end %判断误差是否达到精度4. 计算结果： k,u=kuai_GaussSeidel（n） 1.3280k = 93u = contourf （u, u图二 块Gauss-Sediel迭代法三 （预条件）共轭斜量法求解线性方程组Au=f。1.基本原理（1）预条件共轭斜量法原理（2）预优矩阵的选取2
