matlab Jacobi迭代法Gauss-seidel和SOR迭代

Jacobi与Gauss-Seidel迭代的比较及r算法的MATLAB实现胡志成【摘要】探讨比较了Jacobi迭代和Gauss-Seidel迭代,尤其是它们的可并行性.给出了它们在MATLAB中的一个实现.通过算例说明,算法的并行化实现可以大幅度地提升Jacobi迭代的效率.【期刊名称】《高师理科学刊》【年(卷),期】2018(038)003【总页数】4页(P59-61,84)【关键词】Jacobi迭代;Gauss-Seidel迭代;并行性【作者】胡志成【作者单位】南京航空航天大学理学院,江苏南京 210016【正文语种】中文【中图分类】O241.6;G642.0Jacobi迭代和Gauss-Seidel迭代是解线性方程组的2个经典的迭代法，简单又易于实现，是大学计算方法相关课程必讲的知识点．然而教材[1-3]往往只注重理论而轻视算法的程序实现．就这2个迭代法而言，表现为着重讨论迭代法的收敛性和收敛速度，而甚少给出或者仅给出一个粗糙低效的算法实现．至于Jacobi迭代的并行性，探讨的更是少之又少．由于课时的限制，在工科计算方法教学中更是如此．这一方面是由于一般程序语言关于算法的并行化实现并不那么容易，另一方面则是因为在课堂上短时间内举例说明算法并行化实现带来的效率提升较为困难．然而在课堂上适当地给出算法的实现和具体算例，显然有助于学生对算法本身的理解．本文首先回顾了Jacobi迭代和Gauss-Seidel迭代的基本格式和传统比较，然后探讨了它们的可并行性，并给出了在MATLAB中的一个实现．文中给出的算例不仅可以在课堂的短时间内展示结果，也能说明算法的并行化实现可以大幅度地提升Jacobi迭代的效率．设是阶非奇异方阵，为维列向量．若的主对角元，，则线性方程组可以用Jacobi 迭代和Gauss-Seidel迭代求解．其迭代格式分别为若再将表示成，其中：；；，则Jacobi迭代可用矩阵表示为相应地，Gauss-Seidel迭代的矩阵形式为关于Jacobi迭代和Gauss-Seidel迭代的分析和比较，在众多的文献中都有提及[1，4-6]，主要归结为：（1）收敛性．当且仅当各自的迭代矩阵（分别是和）的谱半径小于1时，该迭代法收敛．若是严格对角占优阵，则2个迭代法都收敛．但一般地， Jacobi迭代收敛不能保证同一问题的Gauss-Seidel迭代也收敛，反之亦然．（2）收敛速度．当Gauss-Seidel迭代和Jacobi迭代两者都收敛时，多数情况下Gauss-Seidel迭代比Jacobi迭代收敛得更快．对于一些模型问题，Gauss-Seidel迭代的收敛速度甚至是Jacobi迭代的2倍以上[5]58．也就是说，从同一初值出发，为得到满足精度要求的解，Gauss-Seidel迭代所需的步数可能比Jacobi迭代所需步数的一半还要少．（3）每步计算量．由于矩阵的求逆并不那么容易，具体的算法实现多采用分量形式，即式（1）和式（2）．比较式（1）和式（2）可知，2个迭代法每步的计算量相等．（4）存储量．Gauss-Seidel迭代只需1个向量来存储或的分量，而Jacobi迭代需要2个向量来分别存储和．比较可知，在传统的算法实现中，Gauss-Seidel迭代比Jacobi迭代往往更具优势．随着多核处理器的迅猛发展，并行计算渐次成为提升计算效率的重要技术手段，因而一个算法的可并行性也成为了衡量其优劣的重要指标[7-9]．进一步比较迭代法（1）和（2），易知Jacobi迭代中，各分量的计算是相互独立的，因而具有高度的并行性；而Gauss-Seidel迭代中，各分量是依次顺序计算的，因而本质上是一个串行算法．这就意味着，利用算法的可并行性，如果实现得当，即使收敛速度较慢，Jacobi迭代仍然有可能在更短的时间内得到满足精度要求的数值解．为了加深学生对此的理解，并让学生体会到算法实现，特别是并行化实现的重要性，有必要给出适当的实例和程序实现．MATLAB是一种面向科学与工程计算的数学软件，因其强大的数值计算和图形绘制能力、广泛的工具包支持以及使用的便利性，深受教学和科研工作者的喜爱，是工科大学生需要掌握的重要软件之一[10]．借助其灵活的矩阵向量运算，给出Jacobi迭代（1）在MATLAB环境中的一个实现（见算法1）和Gauss-Seidel迭代（2）的一个MATLAB实现（见算法2）．可以看到，算法1完全避免了Jacobi迭代各分量之间的for循环，而算法2则无法避免Gauss-Seidel迭代各分量之间的for循环，这是导致2个算法执行效率巨大差异的直接原因．算法1（Jacobi 迭代）：function [x，k，res] = Jacobi（A，b，x0，tol，kmax ）k = 0； res = 1；while（k <= kmax）& res > tol）x =（b - A * x0 + diag（A）．* x0 ）．/ diag（A）；k = k+1；res = max（abs（b - A * x ））；x0 = x；end算法2（ Gauss-Seidel 迭代）：function [x，k，res] = GaussSeidel（A，b，x0，tol，kmax ）k = 0；res = 1；x = x0；while （k <= kmax）&（res > tol）for i = 1：length（x）x（i）=（b（i）- A（i，1：end）* x + A（i，i）* x（i））/ A（i，i）；endk = k+1；res = max（abs（b - A * x））；end以零向量为初值，以残量的无穷范数为收敛准则，考察比较这2个算法求解线性方程组的数值结果，包括迭代步数、每步耗时和总耗时．为了使测试算例具有一定的普遍意义，且能在课堂上较短的时间内展示比较结果．用算法3来随机生成阶线性方程组的系数矩阵．算法3（随机生成一个严格对角占优矩阵）：function A = randSDDMatrix（n）nzmax = randi（[n*5，round（n*10）]，1）；i = [randi（[1，n]，nzmax，1）]；j = [randi（[1，n]，nzmax，1）]；s = rand（nzmax，1）；A = sparse（i，j，s，n，n）；x = randi（[1，1e1]，n，1）+ sum（ abs（A），2）；ii = 1：n；A = sparse（ii，ii，x，n，n）+ A；该矩阵是稀疏的，并且是严格对角占优的，可以确保阶数较大时，2个迭代法仍能快速收敛．当然，右端项也是随机生成的．对每一个阶数，都测试100组和，并取其平均数据（见表1）．由表１可以看出，虽然Jacobi迭代（1）和Gauss-Seidel迭代（2）的每步计算量相等，但是算法1的每步计算效率显著优于算法2．这最终导致：虽然Jacobi 迭代的总迭代步数是Gauss-Seidel迭代的２倍以上，但是总耗时却远少于Gauss-Seidel迭代．而且随着矩阵阶数增加，在本文的算法实现下，这个优势将进一步扩大．本文比较分析了Jacobi迭代和Gauss-Seidel迭代，强调说明Jacobi迭代是一个高度并行的算法，而Gauss-Seidel迭代则是一个串行算法．因而如果算法实现得当，即使收敛速度更慢的时候，Jacobi迭代仍然可能比Gauss-Seidel迭代更快地给出数值解．在教学实践中，课堂上展示这些MATLAB实现和数值算例，明显地加深了学生对这2个迭代法的理解，也使学生意识到了算法实现的重要性，获得了较好的教学效果．【相关文献】[1] 倪勤，王正盛，刘皞．数值计算方法[M]．北京：高等教育出版社，2012[2] 李庆杨，王能超，易大义．数值分析[M]．5版．北京：清华大学出版社，2008[3] 李有法，李晓勤．数值计算方法[M]．2版．北京：高等教育出版社，2005[4] 杜衡吉，徐昆良．Jacobi和Gauss-Seidel迭代法求解线性方程组的分析及应用[J]．曲靖师范学院学报，2011，30（3）：46-50[5] 白红梅．Jacobi迭代法与Gauss-Seidel迭代法的收敛性比较分析[J]．呼伦贝尔学院学报，2009，17（6）：55-58[6] 包丽君．“数值分析”中实践教学的探讨[J]．宁波广播电视大学学报，2008，6（2）：91-93[7] 徐树方，高立，张平文．数值线性代数[M]．北京：北京大学出版社，2000[8] 杨学军．并行计算六十年[J]．计算机工程与科学，2012，34（8）：1-10[9] 韩建勋．并行计算在矩阵运算中的应用[J]．信息与电脑：理论版，2017（5）：93-96[10] 韦慧，蒋利华．基于MATLAB的大学工科《计算方法》课程教学方法探讨[J]．内江科技，2017（4）：72-73。

MATLAB块雅克比迭代法求线性方程组Ax=b的解块高斯-赛德尔迭代法求线性方程组Ax=b的解function [x,N]= BJ(A,b,x0,d,eps,M) %块雅克比迭代法求线性方程组Ax=b的解if nargin==4eps= 1.0e-6;M = 10000;elseif nargin<4errorreturnelseif nargin ==5M = 10000; %参数的默认值endNS = size(A);n = NS(1,1);if(sum(d) ~= n)disp('分块错误！');return;endbnum = length(d);bs = ones(bnum,1);for i=1:(bnum-1)bs(i+1,1)=sum(d(1:i))+1;%获得对角线上每个分块矩阵元素索引的起始值endDB = zeros(n,n);for i=1:bnumendb = bs(i,1)+d(i,1)-1;DB(bs(i,1):endb,bs(i,1):endb)=A(bs(i,1):endb,bs(i,1):endb );%求A的对角分块矩阵endfor i=1:bnumendb = bs(i,1)+d(i,1)-1;invDB(bs(i,1):endb,bs(i,1):endb)=inv(DB(bs(i,1):endb,bs(i ,1): endb));%求A的对角分块矩阵的逆矩阵endN = 0;tol = 1;while tol>=epsx = invDB*(DB-A)*x0+invDB*b; %由于LB+DB=DB-AN = N+1; %迭代步数tol = norm(x-x0); %前后两步迭代结果的误差x0 = x;if(N>=M)disp('Warning: 迭代次数太多，可能不收敛！');return;endendfunction [x,N]= BGS(A,b,x0,d,eps,M) %块高斯-赛德尔迭代法求线性方程组Ax=b的解if nargin==4eps= 1.0e-6;M = 10000;elseif nargin<4errorreturnelseif nargin ==5M = 10000;endNS = size(A);n = NS(1,1);bnum = length(d);bs = ones(bnum,1);for i=1:(bnum-1)bs(i+1,1)=sum(d(1:i))+1;%获得对角线上每个分块矩阵元素索引的起始值endDB = zeros(n,n);for i=1:bnumendb = bs(i,1)+d(i,1)-1;DB(bs(i,1):endb,bs(i,1):endb)=A(bs(i,1):endb,bs(i,1):endb ); %求A的对角分块矩阵endLB = -tril(A-DB); %求A的下三角分块阵UB = -triu(A-DB); %求A的上三角分块阵N = 0;tol = 1;while tol>=epsinvDL = inv(DB-LB);x = invDL*UB*x0+invDL*b; %块迭代公式N = N+1;tol = norm(x-x0);x0 = x;if(N>=M)disp('Warning: 迭代次数太多，可能不收敛！');return;endend类别：matlab 编程 | | 添加到搜藏 | 分享到i贴吧 | 浏览(168) | 评论 (0)上一篇：MATLAB 共轭梯度法求线性方程组A...。

cout<<"发散"<<endl; break; } for(int i=1;i<=N;i++) { float sum1=0; float sum2=0; for( int j=i+1;j<=N;j++) {sum1+=xk[j]*a[i][j];} for( j=1;j<=i-1;j++) {sum2+=a[i][j]*xk[j];} x[i]=(a[i][N+1]-sum1-sum2)/a[i][i]; if(fabs(x[i]-xk[i])<e)
1数值计算方法实验报告五班级地信10801序号姓名实验题目jacobi迭代法和gaussseidel迭代法实验学时
数值计算方法实验报告（五）
班级： 地信 10801
序号：
一、实验题目：jacobi 迭代法和 Gauss-Seidel 迭代法
姓名：
二、实验学时: 2 学时
Hale Waihona Puke 三、实验目的和要求：1．掌握迭代法的基础原理。
times++; Nx=0; if(times>=M){cout<<"发散"<<endl; for(int i=1;i<=N;i++)
break;}
1
{ float sum=0; for(int j=1;j<=N;j++) if(i!=j)sum+=xk[j]*a[i][j]; x[i]=(a[i][N+1]-sum)/a[i][i]; if(fabs(x[i]-xk[i])<e)Nx++;

Jacobi 迭代法与Gauss-Seidel迭代法算法比较目录1 引言 (1)1.1Jacobi迭代法 (2)1.2Gauss-Seidel迭代法 (2)1.3逐次超松弛（SOR）迭代法 (3)2算法分析 (3)3 结论 (5)4 附录程序 (5)参考文献 (8)Jacobi 迭代法与Gauss-Seidel 迭代法比较1 引言解线性方程组的方法分为直接法和迭代法，直接法是在没有舍入误差的假设下，能在预定的运算次数内求得精确解，而迭代法是构造一定的递推格式，产生逼近精确值的序列。

这两种方法各有优缺点，直接法普遍适用，但要求计算机有较大的存储量，迭代法要求的存储量较小，但必须在收敛性得以保证的情况下才能使用。

对于高阶方程组，如一些偏微分方程数值求解中出现的方程组，采用直接法计算代价比较高，迭代法则简单又实用，所以比较受工程人员青睐。

迭代法求解方程组就是构造一个无限的向量序列，使它的极限是方程组的解向量。

即使计算机过程是精确的，迭代法也不能通过有限次算术运算求得方程组的精确解，而只能逐步逼近它。

因此迭代法存在收敛性与精度控制的问题。

迭代法是常用于求解大型稀疏线性方程组（系数矩阵阶数较高且0元素较多），特别是某些偏微分方程离散化后得到的大型稀疏方程组的重要方法。

设n 元线性微分方程组b Ax = (1)的系数矩阵A 非奇异，右端向量0≠b ,因而方程组有唯一的非零解向量。

而对于这种线性方程组的近似解，前辈们发展研究了许多种有效的方法，有Jacobi 迭代法、Gauss —Seidel 迭代法，逐次超松弛迭代法（SOR 法），这几种迭代方法均属一阶线性定常迭代法，即若系数矩阵A 分解成两个矩阵N 和P 的差，即P N A -=；其中N 为可逆矩阵，线性方程组(1)化为：b x P N =-)（b Px Nx +=b N Px N x 11--+=可得到迭代方法的一般公式：d Gx xk k +=+))1(（ （2）其中：P N G 1-=，b N d 1-=，对任取一向量)0(x 作为方程组的初始近似解，按递推公式产生一个向量序列)1(x ，)2(x ，...，）k x(，...，当k 足够大时，此序列就可以作为线性方程组的近似解。