解线性方程组的几种迭代算法

解线性方程组的几种迭代算法内容摘要:本文首先总结了分裂法解线性方程组的一些迭代算法,在此基础上分别通过改变系数矩阵A的分裂形式和对SSOR算法的改进提出了两种新的算法,并证明了这两种算法的收敛性.与其它方法相比,通过改变系数矩阵A的分裂形式得到的新算法具有更好的收敛性,改进的SSOR算法有了更快的收敛速度.最后通过数值实例验证了这两种算法在有些情况下确实可以更有效的解决问题.关键词:线性方程组迭代法算法收敛速度Several kinds of solving linear equationsiterative algorithmAbstract:In this paper, we firstly summarize some Iterative algorithms of Anti-secession law solution of linear equations. Based on these, two new algorithms are put forward by changing the fission form of coefficient matrix A and improving the algorithm of SSOR, and the convergence of the two algorithms is demonstrated. Compared with other methods, the new algorithm acquired by changing the fission form of coefficient matrix A is possessed of a better convergence. And the improved SSOR algorithm has a faster convergence speed. Finally, some numerical examples verify that the two algorithms can solve problems more effectively in some cases.Key words:Linear equations Iteration method algorithm Convergence speed目录1.引言 (1)2.迭代法原理 (1)3.基本迭代法 (2)3.1 Jacobi迭代 (2)3.2 Gauss-Seidel迭代法 (3)3.3 SOR算法 (3)3.4 SSOR算法 (4)3.5 收敛性分析 (4)4.几种新的迭代算法 (5)4.1 基于矩阵分裂形式的新迭代算法 (5)4.2 加权－对称超松弛迭代法 (7)5.算法的不足与改进方法 (9)6.数值实例 (9)6.1渐进收敛速度 (9)6.2几种迭代方法的比较 (10)附录 ............................................... 错误!未定义书签。

参考文献 . (11)解线性方程组的几种迭代算法1.引言在工程技术、自然科学和社会科学中的许多问题最终都可归结为解线性方程组, 因此线性方程组的求解对于解决实际问题是极其重要的.线性方程组的解法有很多种,主要的方法有直接法和迭代法.迭代法就是用某种极限过程去逼近线性方程组精确解的方法,该方法具有对计算机的存贮单元需求少,程序计算简单,原始系数矩阵在计算过程中不变等优点,是求解大型稀疏矩阵方程组的重要方法.目前,人们已经得到了一些较为成熟的线性方程组的迭代解法,从某种意义上讲它们都可归结为分裂法.但在解决具体问题时我们仍面临着许多问题,如:怎样设计出满足要求的求解算法；如何分析、区别算法的好坏；可否改进现有的算法使其更有效；求解所给问题最好可能的算法会是什么,等等.针对这些问题,很多人都做过了大量的研究.文献[2]对迭代法的原理及一些常用的迭代算法进行了研究.文献[1],[3],[4],[5]给出了一些基本的迭代算法并证明了其收敛性.文献[6],[9],[10],[13],[16]研究了一些特殊方程组的迭代解法.文献[7],[8],[12],[14],[15]都是针对不同的问题对超松弛迭代算法进行了改进.文献[11]主要讨论了迭代法解线性方程组的MATLAB 实现.本人在求解线性方程组的问题时,通过对现有迭代算法的改进得到了两种新的算法.本文对这两种算法的收敛性进行了证明,并通过数值实例验证了其在解决某些问题时具有的优势.2.迭代法原理设线性代数方程组为AX b = (1)常常将系数矩阵A 分裂成两个矩阵M 和N 之差,即A M N =- (2)且用迭代(1)1()1k k X M NX M b +--=+ (3)来解线性方程组(1). 将(3)式表示为(1)(),0,1,k k X BX f k +=+= (4)其中,11,B M N f M b --==,称此迭代方法为分裂法,而将B 称为迭代格式(4)的迭代矩阵.然而迭代法需要解决的首要任务是迭代格式是否收敛的问题,任取初始向量()(0)(0)(0)(0)12,,,Tn X x x x =代入(4)中,计算可得迭代序列(1)(2)(1),,.k X X X +若迭代序列{}(1)k X +收敛,设{}()K X 的极限为*X ,对迭代式(4)两边取极限可得:()(1)()lim lim k k k k X BX f +→∞→∞=+即**X BX f =+,*X 是方程组(1)的解,此时称迭代法收敛,否则称迭代法发散.我们有如下的结果:定理2.1[1] 迭代格式(4)收敛的充分必要条件是迭代矩阵B 的谱半径()1B ρ<,而且()B ρ越小,收敛越快.定理2.2[1] 若p B 为矩阵B 的某范数,则总有()pBB ρ≥.对于矩阵A 的分裂应该说是有很多形式,但并不是所有分裂形式产生的迭代格式都有意义.于是我们有正规分裂的概念:定义2.1[2]对于实方阵A ,若矩阵M 和N 满足A M N =-,且10,0M N N -≥≥则称A M N =-是A 的一个正规分裂.那么正规分裂与其他分裂形式相比到底有什么优势呢?我们有如下定理:定理2.3[2]若A M N =-为A 的正规分裂,且10A -≥,则()()()1111A N M N A N ρρρ---=+从而()11M N ρ-<,此时相应的迭代格式(3)必收敛.如果针对矩阵A 给出两种正规分裂,如何来衡量它的好坏呢?定理2.4[2] 若矩阵有两个正规分裂,设1122,A M N A M N =-=-且10A -≥,211120,0,N N N N N ≥≥≠≠,则有()()11112101M N M N ρρ--≤≤<.3.基本迭代法下面给出常见的几种基本迭代格式.将()ij A a =分裂为A D L U =-- (5)其中,11121212221,10000,,00n n n n n nn a a a a a a D L U a a a -⎛⎫⎛⎫⎛⎫⎪⎪⎪ ⎪⎪⎪==-=- ⎪ ⎪⎪ ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭ 3.1 Jacobi 迭代取,.M D N L U ==+ (6)则(4)式中迭代矩阵B 和右端向量f 分别为()11,.B D L U f D b --=+=迭代格式(4)的分量形式为(1)()1,1,2,.0,1,.k k ii ij j j i ii xb a x i n k a +≠⎛⎫=-== ⎪⎝⎭∑称它为Jacobi 迭代法,该迭代法具有和(1)k x +中分量的计算次序无关,容易并行计算等优点.3.2 Gauss-Seidel 迭代法取,,M D L N U =-= (7)此时(4)式中迭代矩阵B 和右端向量f 分别为()()11,.B D L U f D L b --=-=-迭代格式(4)的分量形式为(1)(1)()111,1,2,.0,1,.k k k i i ij j ij j j i j i ii x b a x a x i n k a ++≤-≥+⎛⎫=--== ⎪⎝⎭∑∑称它为Gauss-Seidel 迭代法.和Jacobi 迭代法相比Gauss-Seidel 迭代法使用了最新已经计算的分量. 3.3 SOR 算法取11,,M D L N M A D U ωωω-=-=-=+ (8)此时(4)式中迭代矩阵B 和右端向量f 分别为()()()()111,.B D L D U f D L b ωωωωω--=--+=-迭代格式(4)的分量形式为()1(1)()(1)()111,1,2,.0,1,i nk k k k iii ij jij j j j i ii xxb a x a x i n k a ωω-++==+⎛⎫=-+--== ⎪⎝⎭∑∑其中,ω称为松弛因子(总假定是实数).迭代格式的矩阵形式为:()()()11(1)()1k k X D L D U X D L b ωωωωω--+=--++-⎡⎤⎣⎦称它为对应于松弛因子ω的逐次超松弛迭代法(Successive Over Relaxation,简称SOR).它可以看成Gauss-Seidel 迭代法与原向量的组合,但使用了最新已经计算的分量.也就是把(1)k i x +取为Gauss-Seidel 迭代法中()k i x 与(1)k i x +的某个平均值,即:()()(1)()(1)()(1)()1k k k ii i k k k iiix x x xxxωωω+++=-+=+-当1ω=时,它就是Gauss-Seidel 迭代法.因此希望选取合适的ω使得它比Gauss-Seidel 迭代法具有更快的收敛速度. 3.4 SSOR 算法SOR 迭代格式还可以写为()()(1)()1,0,1,k k D L X D U X b k ωωωω+-=-++=⎡⎤⎣⎦将上式L 和U 的位置互换就得到一个新的迭代格式,具体表示为:()()()()1()()21()(1)211,0,1,k k k k D L X D U X b D L X D U X b k ωωωωωωωω+++⎧-=-++⎡⎤⎪⎣⎦⎨⎪-=-++=⎡⎤⎣⎦⎩若消去1()2k X+,就得到迭代格式(4),其中()()()()()()()1111111,2B D U D D L D L D D U f D U D D L ωωωωωωωωωω-----=---+-+⎡⎤⎡⎤⎣⎦⎣⎦=---称为对称逐次超松弛迭代法(Symmetric Successive Over Relaxation ,简称SSOR). 对一类椭圆微分方程离散后得到的线性方程组,Young[3]给出了最佳松弛因子,即为ω=其中B 是Jacobi 迭代法中的迭代矩阵.在实际问题中最佳松弛因子是很难计算的,但一般都在(0,2)之间. 3.5收敛性分析定义3.1[3] 若实矩阵()ij nnA a =,满足()1,,1,2,nij ii j j ia a i n =≠<=∑或1,,nij jj i i ja a =≠<∑()1,2,j n =则称A 为严格对角占优矩阵；若满足()1,,1,2,nij ii j j ia a i n =≠≤=∑或1,,(1,2,)nij jj i i ja a j n =≠≤=∑且上述不等式至少有一个严格成立,则称A 为弱严格对角占优矩阵.定义3.2[3] 设A 为n 阶矩阵,2n ≥.若存在n 阶排列矩阵P ,使得1112220T AA P AP A ⎛⎫= ⎪⎝⎭其中11A 为r 阶矩阵,1r n ≤<,则称A 是可约的,否则称A 不可约.定理3.1[3] 若A 为n 阶严格对角占优矩阵或不可约的弱严格对角占优矩阵,则对任意的初始值(0)x ,(1)式中的Jacobi 迭代法、Gauss-Seidel 迭代法和关于01ω<≤的SOR 迭代法均收敛.定理3.2 [4] 对所有ω均成立不等式()1SOR B ρω≥-当ω是实数时,SOR 方法收敛的一个必要条件是02ω<<.定理3.3[4] 如果系数矩阵A 是Hermite 矩阵,则SOR 方法收敛的充要条件是:A 正定和02ω<<.4.几种新的迭代算法4.1 基于矩阵分裂形式的新迭代算法上述几种基本迭代方法都是通过对线性方程组(1)的系数矩阵A 进行分裂得到的,不同之处在于A 分裂成A M N =-的形式时,M 和N 的取值不同.由(3),(4)可知此种格式下的迭代矩阵11B M N E M A --==-,所以当M 是A 一个很好的近似时,p B 就会很小.再由定理2.1和定理2.2可知此时得到的迭代法收敛速度也更快.另一方面,我们构造的迭代格式为(1)(),0,1k k Mx Nx b k +=+=,由于每一次迭代都要解一个方程组,所以我们也要求非退化的矩阵M 的形式比较简单,如对角矩阵、下三角矩阵等.比较Jacobi 迭代、Gauss-Seidel 迭代及SOR 算法中M 的取值,我们可以看到M 的形式都为对角矩阵或下三角矩阵,并且随着M 越近似等于A ,所得到的迭代方法的收敛速度也越快.满足上述两个条件,对M 取不同于这三种算法中M 取值的形式,我们可以得到下述新的算法.4.1.1 算法的建立 取(),1.01M D L N L U ωωω=-=-+≤≤, (9)作为A M N =-的一个新分裂,在此分裂的基础上可得到一个新的迭代公式(1)(),0,1,k k x B x f k ωω+=+= (10)其中,()()()11B D L L U ωωω-=--+,()1,01f D L b ωωω-=-≤≤.显然,当0ω=时就是Jacobi 迭代法；当1ω=时就是Gauss-Seidel 迭代法.4.1.2 收敛性分析引理1[3] A 可约的充要条件为存在非空子集{}1,2,J n ⊂,使得0,,.ij a i J j J =∈∉引理2[3] 若A 为n 阶严格对角占优矩阵或不可约的弱严格对角占优矩阵,则A 是非退化的.定理 4.1 若A 为n 阶严格对角占优矩阵或不可约的弱严格对角占优矩阵,则对任意的初始值(0)x ,当01ω≤≤时,(10)式表示的迭代格式收敛.证明:(10)式的迭代矩阵为()()()11B D L L U ωωω-=--+我们有()()()()()()()()()()()11det det 1det det 1E B E D L L U D L D L L U ωλλωωωλωω---=---+=----+其中λ为B ω的特征值.若A D L U =--为严格对角占优矩阵,则D L ω-也为严格对角占优矩阵；若A D L U =--为不可约的弱对角占优矩阵,则D L ω-也为不可约的弱对角占优矩阵.由引理2,我们有()()1det 0D L ω--≠所以()()()()()det 0det 10E B D L L U ωλλωω-=⇔---+=令()()()()()()()()11121312122232123111111n n n n n nn H D L L U D L Ua a a a a a a a aa a a λωωλλωωλλωωλλωωλωωλωωλ=---+=--+-⎛⎫⎪-+ ⎪= ⎪ ⎪ ⎪-+-+-+⎝⎭下面证明1λ<,用反证法.假设1λ≥,由01ω≤≤可知A 和H 的对应元素一定同时为0或非0.由引理1知,A 不可约推出H 也不可约.下面比较λ和1λωω-+的大小. ①当1λ≥时,()()()()11111110λλωωλωλλωλλω--+=--+=---=--≥②当1λ≤-时,110λ--≤≤,()()()()()1111111110λλωωλλωωλλλωλλλωλ---+=--+-+=---+=---+≥由①、②知:当1λ≥,01ω≤≤时,1λλωω≥-+.所以有:A 为弱对角占优矩阵时,H 也是弱对角占优矩阵；A 为严格对角占优矩阵时,H 也是严格对角占优矩阵.再由引理2,可以得到H 为非退化的矩阵.从而,()det 0H ≠.进而,()det 0E B ωλ-≠.这与λ是B ω的特征值矛盾,所以1λ<.进一步可得到()1B ωρ<,我们就证明了(10)式表示的迭代格式收敛.■4.2 加权－对称超松弛迭代法4.2.1 算法的建立考虑SOR 算法和4.1给出的新迭代算法,其实质都是给矩阵A 分裂以后的因子一个权重后得到的,目的是使得迭代格式(4)的收敛速度尽可能快.按照这种思想,如果在SOR 和SSOR 迭代中引入一个参数θ,就可以得到迭代格式()()()(1)()()1,[0,1]k k k SOR SOR SSOR SSOR X B X f B X f θθθ+=++-+∈ (11) 其中,SOR B ,SSOR B 分别为SOR 迭代和SSOR 迭代算法的迭代矩阵.此时的迭代矩阵()1SOR SSOR B B B θθ=+-.显然,当0θ=时,该算法即为SSOR 迭代；当1θ=时,即为SOR 迭代.若迭代格式(11)收敛到某个向量*X ,则有()()()***1,[0,1]SOR SOR SSOR SSOR X B X f B X f θθθ=++-+∈所以不管θ在[0,1]内如何取值,当迭代格式(11)收敛时,它必定是原来方程X BX f =+的解.选取适当的B ,根据(4)式的定义,这样的解也是方程(1)的解.显然我们可以选择一个较好的θ,使得迭代格式(11)收敛的尽可能快.已知SOR 迭代的核心部分为:()1(1)()(1)()111:1:1i nk k k k iii ij jij j j j i ii for i nx xb a x a x a endωω-++==+=⎛⎫=-+--⎪⎝⎭∑∑ SSOR 迭代的核心部分为:()()1(1)()(1)()111(1)()(1)()111:1:1:1:11i nk k k k iii ij jij j j j i ii i nk k k k iii ij jij j j j i ii for i nxxb a x a x a endfor i n x xb a x a x a endωωωω-++==+-++==+=⎛⎫=-+--⎪⎝⎭=-⎛⎫=-+--⎪⎝⎭∑∑∑∑将SOR 迭代和SSOR 迭代生成的向量进行加权平均,可得到改进的SSOR 迭代算法,其核心部分为:()()()1(1)()(1)()111(1)()(1)()11(1)(1)(1)1:1:1:1:111:1:1i nk k k k iii ij jij j j j i ii i nk k k k i ii ij jij j j j i ii k k k iii for i nx xb a x a x a endfor i n xxb a x a x a endfor i nx x x endωωωωθθ-++==+-++==++++=⎛⎫=-+--⎪⎝⎭=-⎛⎫=-+--⎪⎝⎭==+-∑∑∑∑ 其中,θ为权重,可以取[0,1]之间的任意实数,在此称之为加权－超松弛迭代法.4.2.2 收敛性分析定理4.2 加权－超松弛迭代算法,满足迭代关系式(4),其迭代矩阵()1SOR SSOR B B B θθ=+-.其中,()()11SOR B D L D U ωωω-=--+⎡⎤⎣⎦()()()()11111SSOR B D U D D L D L D D U ωωωωωω---=---+-+⎡⎤⎡⎤⎣⎦⎣⎦ 且(11)对任意的初始向量都收敛的充要条件为()1B ρ<.证明:由4.2.1中的分析知该定理结论显然成立. ■ 5.算法的不足与改进方法上述两种算法本质上都是添加了一个加权因子,不同之处在于一个是对分裂后的矩阵添加的,而另一个是对迭代矩阵添加的.增加了加权因子之后,如何对其取值才最合理就变成了一个亟需解决的问题.像SOR 算法中很难确定ω的值一样,我们很难找到一种通用的方法来确定所增加的加权因子的值.在下面的数值实例中,为了验证算法的优越性,我们可以让加权因子ω遍取[0,1]中的某些值,通过MATLAB 编程搜索出迭代次数最少时ω的取值.但是在解决实际问题时如果仍这样做,不但不能减少迭代的总次数,反而是增加了计算机的工作量和解决问题所需要的时间.下面提供一种解决实际问题时加权因子ω值的确定方法:(1)给定某个迭代次数n ,搜索出迭代n 次后解的精度最小的那个ω的值；（2）用（1）中ω的值作为加权因子的值求解线性方程组的指定精度的解.6.数值实例6.1渐进收敛速度为方便比较和表述用不同的迭代方法解线性方程组时的收敛速度,我们引入下述定义.定义6.1[5]称()1ln k k k R B B=-为迭代法()(1)k k X BX f -=+的平均收敛率. 此处定义的是平均压缩率的对数值(再取负号),它依赖于所选择的范数和迭代次数,使用起来非常不便.由于()1lim k k k B B ρ→∞=则再给出下面的定义:定义6.2[5]称()()ln R B B ρ=-为迭代法的渐进收敛速度,简称迭代法的收敛速度.6.2几种迭代方法的比较下面我们用不同的迭代方法求解n 元线性方程组Ax b =,其中4131412,.1412143A b -⎛⎫⎛⎫ ⎪ ⎪-- ⎪ ⎪ ⎪ ⎪== ⎪ ⎪-- ⎪ ⎪ ⎪ ⎪-⎝⎭⎝⎭方程的精确解为()*1,1,,1Tx =.取10n =,初始向量(0)x 为零向量.下表给出了几种迭代方法达到不同精度时所需的迭代次数.(程序见附录)程组时,4.2中提供的方法收敛速度最快.根据6.2的定义,我们算出精度为8110-⨯时各种迭代算法的渐进收敛速度,见表2. (程序见附录)表2 精度为8110-⨯时各种迭代算法的渐进收敛速度 结合表1与表2,我们可以看到:当迭代结果的精度达到8110-⨯时,表1显示的Jacobi 迭代法、Gauss-Seidel 迭代法、SOR 算法、SSOR 算法和对称－超松弛迭代算法所需的迭代次数逐渐减少,4.2中给出的方法是收敛最快的,这与表2中给出的其渐进收敛速度逐渐变大的结果一致.我们还可以看到,用4.1中的算法进行计算所需的迭代次数较SSOR 方法的迭代次数多,但由定义6.2计算出的其渐进收敛速度却较小,这说明用该定义意义下的范数去度量4.1中的算法的收敛速度,并不是一个好的度量方式.考虑方程组123101*********x x x ⎛⎫⎛⎫⎛⎫ ⎪⎪ ⎪=- ⎪⎪ ⎪ ⎪⎪ ⎪--⎝⎭⎝⎭⎝⎭-1其精确解为()2,5,3T-.求出各种迭代方法的迭代矩阵B 对应的谱半径()B ρ,可以看到只有Jacobi 迭代法和4.1中的算法是收敛的,其他方法均发散.下面我们比较Jacobi 迭代法和4.1中的算法在达到指定精度时的收敛次数,见表3.(程序见附录)表3 两种迭代法达到指定精度所需迭代次数的比较参考文献[1] 张韵华,奚梅成,陈效群,数值计算方法与算法[M].北京:科学出版社,2009: 130-140.[2] 邵新慧,大型线性方程组的迭代解法[J ].东北大学博士学论文,2009.[3] 李大明,数值线性代数[M].北京:清华大学出版社,2010:177-192.[4] 余德浩,汤华中.微分方程数值解法[M].北京:科学出版社,2004:320-328.[5] 薛毅,数值分析与实验[M].北京:北京工业大学出版社,2005:88-93.[6] 马云,数值分析中的迭代法解线性方程组[J ].考试周刊,2010:71-72.[7] 段班祥,吴教育,朱小平,非线性互补问题的改进超松弛迭代算法[J].江西师范大学学报.2009,33(5):617-621.[8] 蒋家羚,王勇,最优超松弛因子的一种确定方法及其在裂纹计算中的应用[J].机械强度,2002,24(1):133-135.[9] 郑亚敏,迭代法解线性方程组的收敛性比较[J ].江西科学.2009.10: 659-661.[10] 谢刚,不定线性方程组的一种迭代算法[J].科学技术与工程.2007.第7卷,第14期.[11] 张步林,线性代数方程组迭代解法的MATLAB 实现[J].成都纺织高等专科学校学报,2008.10:45-47.[12] 李春光,徐成贤,确定SOR最佳松弛因子的一个实用算法[J].计算力学学报.2002.8:299-302.[13]曾闽丽,加权-对称超松弛算法[J].莆田学院院报,2008.4,第15卷,第2期.[14] 周荣富,蔡治亚,关于超松弛迭代法收敛的一个判别准则[J].1991.6.第6卷.第2期.[15] 唐永建,向隆万,最佳双松弛因子迭代法的存在条件[J].上海电力学院学报,第17卷,第1期.[16] 高虹霓,曹泽阳,一种新的非线性方程求根迭代法[J].空军工程大学学报,2002.4.第3卷,第2期.。