QQ(1)(k)?IHk?Q(k?1)
i?1,2,?,n
ti?1n?k?qj?k(k)(k)iju?k?j j?k,k?1,?,nqij(k?1)
?k?j?qij?tiu(5)按要求输出,结束
3．变量说明：
A － 输入的n阶矩阵，同时用于存储上三角阵R； n Q
－ 矩（方）阵A的阶数；
－ 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），对应初等反射阵的ρ,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; %利用向量内积求和 Q(ii,k:n)=Q(ii,k:n)-t2*u'; end end
5．使用示例: (1)A为3阶矩阵： >> A=[1 2 3; 2 3 0; 3 4 5] A =
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.4082 r =
-3.7417 -5.3452 -4.8107 0 0.6547 0.4364 -0.0000 0.0000 3.2660 检验： >> q*r
ans =
1.0000 2.0000 3.0000 2.0000 3.0000 0.0000
3.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] A =
1 2 3 4 2 3 0 1 3 4 5 6 1 6 8 0
>> [q,r]=qrhh(A)
q =
-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 r =
-3.8730 0 0 0 检验： >> q*r
ans =
1.0000 2.0000 3.0000 1.0000
-6.7132 -6.7132 -6.1968 4.4647 6.4805 -1.4783 -0.0000 -3.3070 -3.0178 0.0000 0 -1.8187 2.0000 3.0000 4.0000 3.0000 -0.0000 1.0000 4.0000 5.0000 6.0000 6.0000 8.0000 0
数值求解正方形域上的Poisson方程边值问题
???2u?2u????2?2??f(x,y)?xy,0?x,y?1?y????x??u(0,y)?u(1,y)?u(x,0)?u(x,1)?1
用MATLAB语言编写求解此辺值问题的算法程序，采用下列三种方法，并比较三种方法的计算速度。1、用SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子；2、用块 Gauss-Sediel迭代法求解线性方程组Au=f；3、（预条件）共轭斜量法。 用差分代替微分，对Poisson方程进行离散化，得到五点格式的线性方程组
4ui,j?ui?1,j?ui?1,j?ui,j?1?ui,j?1?hfij2?i,j?Nu1,j?uN?1,j?ui,1?ui,N?1?1，写成矩阵形式Au=f。其中 ?A212??IA???????IA33????I??v2????v3??u??????I????vAN?1,N?1??N??T2 1?i,j?N?1
?b2?b3f??????bN?4?1????14??Aii?????????1??????1??4?v2?(u2,2,u3,2,...,uN,2),v3?(u2,3,u3,3,...,uN,3),......,vN?(u2,N,u3,N,...,uN,N)22TTTTb2?h(f2,2,f3,2,...,fN,2)?(u2,1?u1,2,u3,1,...,uN,1?uN?1,2),b3?h(f2,3,f3,3,...,fN,3)?(u1,3,0,...,fN?1,3),......,bN?h(f2,N,f3,N,...,fN,N)?(u2,N?1?u1,N,u3,N?1,...,uN,N?1?uN?1,N)h?fi,j
2TTTT1N?xiyj,,取N?10,则h?0.1i,j?2,3,...,N
一 用SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子。
1. 基本原理：
Gauss-Seidel迭代法计算简单，但是在实际计算中，其迭代矩阵的谱半径常接近1，因此收敛很慢。为了克服这个缺点，引进一个加速因子（又称松弛因子）对Gauss-Seidel方法进行修正加速。
假设已经计算出第k步迭代的解（i=1，2，···，n），要求下一步迭代的解（i=1，2，···，n）。首先，用Gauss-Seidel迭代格式计算x(k?1)
? x(k?1)??a12x(k)?a13x(k)???a1nx(k)?b12n?1a11a113a11a11?
?(k?1)a23(k)a2n(k)a21(k?1)b2x??x?x???x??13n 2a22a22a22a22??? ?an1(k?1)an2(k?1)ann?1(k?1)bn?(k?1) x??x?x???x?12n?1?naaaannnnnnnn?
然后引入松弛因子，用松弛因子对和作一个线性组合。
xi(k?1)??xi(k?1)?(1??)xi(k)，i=1,2,…,n
将二者合并成为一个统一的计算公式：
(k?1)(k)(k)(k)?a11x1??(?a12x2???a1nxn?b1)?(1??)a11x1?(k?1)(k?1)(k)(k)(k) ??(?a21x1?a23x3???a2nxn?b2)?(1??)a22x2?a22x2? ??(k?1)(k)?ax(k?1)??(?ax(k?1)?ax(k?1)???axn?1?bn)?(1??)annxnnnnn11n22nn?1?
2. 算法
(1)Gauss-Seidel迭代法引入松弛因子w： (k?1)(k?1)(k)ui,j??ui,j?(1??)ui,j,五点格式即为：
du?ui,j?(k?1)
i?2,3,?,n?ui,j(k)?4(bi,j?ai?1,jui?1,j?ai,j?1ui,j?1?ai,j?1ui,j?1?ai?1,jui?1,j?ai,jui,j)(k?1)(k?1)(k)(k)(k)(2)计算步骤：
第一步：给松弛因子赋初值w=1.1～1.8,给场值u和场源b赋初值 第二步：用不同的w进行迭代计算。置error=0;
计算
ui,j(k?1)?ui,j?(k)?4(bi,j?ui?1,j?ui,j?1?ui,j?1?ui?1,j?4ui,j(k?1)(k?1)(k)(k)(k?1))