当前位置:
文档之家› Jacobi G-S SOR迭代法在matlab中的例子
Jacobi G-S SOR迭代法在matlab中的例子
4 / 15
[数值分析报告]
[数学基地班 赵晨晓 2011301000007]
x0=zeros(20,1),运行程序,得到i=38。 Jacobi迭代法绘图程序: function y=jacobian1(A,b,x0) residue=zeros(1000,1); if size(b,2)~=1 b=b'; end if size(x0,2)~=1 x0=x0'; end D=diag(diag(A)); B=D\(D-A); g=D\b; x1=B*x0+g; i=1; while(norm(x1-x0,inf)>1e-10) residue(i)=norm(x1-x0,inf); x0=x1; x1=B*x1+g; i=i+1; end semilogy(residue);
(k)
已经算出,用它代替 x1(k-1) ,其他分量仍用 xi(k-1) 。类似
的,计算 xl(k)时,因 x1(k) , , ,xl-1(k)都已算出,用它们代替 x1(k-1) , , , ,xl-1(k-1) ,其他分量仍用 xk-1 的分量,于是有 Xk=D-1Lxk+D-1Uxk-1+g, k=1,2, , , (6)
10 / 15
[数值分析报告]
[数学基地班 赵晨晓 2011301000007]
end D=diag(diag(A)); L=-tril(A-D); U=-triu(A-D); M=(D-w*L)\((1-w)*D+w*U); g=w*(D-w*L)\b; x1=M*x0+g; i=1; while(norm((x1-x0),2)>1.0e-6) residue(i)=norm(x1-x0,inf); x0=x1; x1=M*x0+g; i=i+1; end semilogy(residue); y=x1; i 由公式 wopt=
3 / 15
[数值分析报告]
[数学基地班 赵晨晓 2011301000007]
if size(b,2)~=1 b=b'; end if size(x0,2)~=1 x0=x0'; end D=diag(diag(A)); B=D\(D-A); g=D\b; x1=B*x0+g; i=1; while(norm((x1-x0),2)>0.0001) x0=x1; x1=B*x1+g; i=i+1; end y=x0; i-1 其中,A是一个20阶矩阵。由于Jacobi迭代法需满足条件A及 2D-A均正定,故可取A=rand(20),令A=A*A’,然后令 A=A+100*eye (20) , 使A满足严格对角占优。 算得cond (A , inf)=3.1223,说明该矩阵条件数较好。令b=ones(20,1),
B <1 或者 B
1
2、 设 B 是 jacobi 迭代的迭代矩阵。 若 则 G-S 迭代收敛。
<1,
3、 若系数矩阵 A 对称, 对角线元素 aii>0 (i=1, 2, , , n) , 则 jacobi 迭代收敛的充分必要条件是 A 和 2D-A 都正定;若 A 正定,则 G-S 迭代收敛。 三、MATLAB 中的矩阵实验 Jacobi 矩阵算法: function y=jacobian(A,b,x0)
[数学基地班 赵晨晓 2011301000007]
可以看出图像弯曲不平,性质没有上一条好。
可以看出, 只需要 6 次迭代, 就能满足要求。 因为 w 较小时, 收敛稳定性高,w 较大时,收敛速度较快。 另一方面,当矩阵条件数没有那么好的时候,例如可令 A=rand(20) ,A=A*A’, cond(A,inf)= 2.8346e+04,此时, 由于 jacobi 迭代条件较为苛刻,jacobi 迭代图像如下所示:
5 / 15
[数值分析报告]
[数学基地班 赵晨晓 2011301000007]
y=x0; i 可 得 到 此 矩 阵 的 收 敛 图 像 :
G-S迭代法算法: function y=GS(A,b,x0) if size(b,2)~=1 b=b'; end if size(x0,2)~=1 x0=x0'; end
9 / 15
[数值分析报告]
[数学基地班 赵晨晓 2011301000007]
L=-tril(A-D); U=-triu(A-D); M=(D-w*L)\((1-w)*D+w*U); g=w*(D-w*L)\b; x1=M*x0+g; i=1; while(norm((x1-x0),2)>1.0e-6) residue(i)=norm(x1-x0,inf); x0=x1; x1=M*x0+g; i=i+1; end y=x1; i SOR 迭代法绘图算法: function y=SOR(A,b,x0,w) residue=zeros(1000,1); if size(b,2)~=1 b=b'; end if size(x0,2)~=1 x0=x0';
13 / 15
[数值分析报告]
[数学基地班 赵晨晓 2011301000007]
这说明此时 jacobi 迭代法不收敛; 而 G-S 迭代法图像如下所示:
14 / 15
[数值分析报告]
[数学基地班 赵晨晓 2011301000007]
这说明它需要迭代 3000 多次才能达到预期目标。 而对于 SOR 迭代法,由于ρ(B)过大,它不收敛。 综上,可以得到结论: (1) 当矩阵收敛时, G-S 迭代法收敛速度总比 Jacobi 迭代法收敛速度要快。而选取恰当 w 后,SOR 迭代法速度比 G-S 迭代法收敛速度 还要快。 (2) 当矩阵条件数较小时,即此时矩阵性质比较 好,三种算法的收敛速度都比较快。 (3) 以上启示我们在进行矩阵迭代之前,最好对 其进行预处理,以使其条件数较小,性质优 良。
6 / 15
[数值分析报告]
[数学基地班 赵晨晓 2011301000007]
D=diag(diag(A)); L=-tril(A-D); U=-triu(A-D); M=(D-L)\U; g=(D-L)\b; x1=M*x0+g; i=1; while(norm((x1-x0),2)>1.0e-6) x0=x1; x1=M*x0+g; i=i+1; end y=x1; i-1 同样的 A,得到 i=8。 C-G 迭代法绘图算法: function y=GS(A,b,x0) residue=zeros(1000,1); if size(b,2)~=1 b=b'; end if size(x0,2)~=1
1 1 1 ( B)* ( B)
利用函数 vrho,代入计算,得到
wopt=0.9,此时收敛图像为:
11 / 15
[数值分析报告]
[数学基地班 赵晨晓 2011301000007]
可以看出这是一条较为光滑且下降较快的曲线。作为对比, 当 w=1.5 时,其图像为:
12 / 15
[数值分析报告]
其中,D 为严格下三角阵,L 为严格上三角阵,那个(1)可 写为 x=Bx+g (4)
其中,B=D-1(L+U),g=D-1b.若给定初始向量 X0=(x1(0) , x2 ( 0 ) , , , , x n( 0 ) )T,并代入(4)的右端,就 可以计算出一个新的向量 x1,即 x1=Bx0+g;
[数值分析报告]
[数学基地班 赵晨晓 2011301000007]
关于 Jacobi G-S SOR 方法的收敛速度比较 摘要: 本论文主要通过判断比较 Jacobi G-S SOR 三种
计算方法的收敛速度,来分析三种算法的优劣。辅助软件为 MATLAB。 随着计算技术的发展,计算机的存储量日益增大,计算速 度也迅速提高,直接法在计算机上可以求解的线性方程组的 规模也越来越大, 但直接法大多需要对系数矩阵 A 进行分解, 因此一般不能保证 A 的稀疏性。而实际应用中,特别是偏微 分方程的数值求解时,常常遇到的就是大型稀疏线性方程的 求解问题。因此寻求能够保持稀疏性的有效算法就成为数值 线性代数中的一个非常重要的研究课题。 一、三种迭代方法内容简介 Jacobi 迭代 考虑非奇异线性代数方程组 Ax=bHale Waihona Puke 令 A=D-L-U (1) (2 )
1 / 15
[数值分析报告]
[数学基地班 赵晨晓 2011301000007]
再把 x1 代入(4)右端,又可得到一个向量 x2;以此类推有 Xk=Bxk-1+g,k=1,2,3, , , (5 )
这被称为 jacobi 迭代格式。B 称为 jacobi 迭代的迭代矩阵,g 称为常数项。 Gauss-Seidel 迭代 假设不按 jacobi 格式,而是在计算 xk 的第一个分量用 xk-1 的各个分量计算。但当计算 xk 的第二个分量 x2(k)时,因 x1
k+1
-xk,则有
xk+1= x+xk。
Xk+1 可以看作是在向量 xk 上加上修正项 x 而得到的。 若
2 / 15
[数值分析报告]
[数学基地班 赵晨晓 2011301000007]
修正项的前面加上一个参数ω,便得到松弛法的计算公式 Xk+1=xk+ω x=(1-ω)xk+ω(D-1Lxk+1+D-1Uxk+D-1b) 其中ω叫做松弛因子。当ω>1 时叫做超松弛;ω<1 时叫 做低松弛; ω=1 即为 G-S 迭代。 我们把超松弛迭代简称为 SOR。
可得到此矩阵的收敛图像:
8 / 15
[数值分析报告]