数值计算-第4章--解线性方程组的迭代法.doc

数值计算_第4章 解线性方程组的迭代法 作者： 日期：30 个人收集整理 勿做商业用途第4章 解线性方程组的迭代法用迭代法求解线性方程组与第4章非线性方程求根的方法相似，对方程组进行等价变换，构造同解方程组(对可构造各种等价方程组，如分解，可逆，则由得到），以此构造迭代关系式 （4。1）任取初始向量，代入迭代式中，经计算得到迭代序列。若迭代序列收敛，设的极限为，对迭代式两边取极限即是方程组的解,此时称迭代法收敛，否则称迭代法发散。我们将看到，不同于非线性方程的迭代方法，解线性方程组的迭代收敛与否完全决定于迭代矩阵的性质，与迭代初始值的选取无关。迭代法的优点是占有存储空间少,程序实现简单,尤其适用于大型稀疏矩阵；不尽人意之处是要面对判断迭代是否收敛和收敛速度的问题.可以证明迭代矩阵的与谱半径是迭代收敛的充分必要条件，其中是矩阵的特征根。事实上，若为方程组的解，则有再由迭代式可得到由线性代数定理，的充分必要条件.因此对迭代法（4。1）的收敛性有以下两个定理成立。定理4.1 迭代法 收敛的充要条件是。定理4.2 迭代法收敛的充要条件是迭代矩阵的谱半径因此,称谱半径小于1的矩阵为收敛矩阵。计算矩阵的谱半径，需要求解矩阵的特征值才能得到，通常这是较为繁重的工作.但是可以通过计算矩阵的范数等方法简化判断收敛的工作。前面已经提到过，若|Ap矩阵的范数，则总有。因此，若，则必为收敛矩阵。计算矩阵的1范数和范数的方法比较简单，其中于是,只要迭代矩阵满足或，就可以判断迭代序列是收敛的。要注意的是,当或时，可以有,因此不能判断迭代序列发散。在计算中当相邻两次的向量误差的某种范数小于给定精度时,则停止迭代计算,视为方程组的近似解（有关范数的详细定义请看3。3节。）4.1雅可比（Jacobi）迭代法4。1.1 雅可比迭代格式雅可比迭代计算元线性方程组（4.2）写成矩阵形式为。若将式( 4。2）中每个方程的留在方程左边,其余各项移到方程右边；方程两边除以则得到下列同解方程组:记,构造迭代形式或 (4。3）迭代计算式（4。3）称为简单迭代或雅可比迭代。任取初始向量 ，由式（4.3）可得到迭代向量序列雅可比迭代矩阵设由，得到等价方程：记不难看出,正是迭代式（4.3）的迭代矩阵,是常数项向量.于是式（4.3)可写成矩阵形式：（4。4)其中：雅可比迭代算法下面描述解线性方程组的雅可比迭代算法，为了简单起见,在算法中假定矩阵满足雅可比迭代要求，即，并设由系数矩阵构造迭代矩阵是收敛的。1定义和输入系数矩阵与常数项向量的元素。2FOR i：=1,2,，n/假定，形成常数项向量 FOR j:=1，2,，n /形成迭代矩阵元素3 / 赋初始值，x1和x2分别表示和4WHILE x1：=x2 x2：=Bx1+g/ FOR u：=1，2，,n/ s:= gu；/FOR v:=1，2,，n s：=s+buvx1v；/ x2u:=s； ENDWHILE5输出方程组的解例4.1 用雅可比方法解下列方程组:解：方程的迭代格式：或雅可比迭代收敛.取初始值,计算结果由表4。1所示。表4.1 计算结果0111 11。51。60.90.2521.252。081.090.483-0。9152.0681。0170。33540.95751.98640.98470。081651。014451。988440.997110.056956-1。007222。002311.00260.0138770.9975432.001971。000490。009687方程组的准确解是4。1.2 雅可比迭代收敛条件对于方程组，构造雅可比迭代格式其中，。当迭代矩阵的谱半径时，迭代收敛，这是收敛的充分必要条件。迭代矩阵的某范数时，迭代收敛。要注意的是范数小于1只是判断迭代矩阵收敛的充分条件，当迭代矩阵的一种范数|B|1，并不能确定迭代矩阵是收敛还是发散。例如，则，但它的特征值是0.9和0.8。是收敛矩阵。当方程组的系数矩阵具有某些性质时，可直接判定由它生成的雅可比迭代矩阵是收敛的.定理4。3 若方程组的系数矩阵，满足下列条件之一，则其雅可比迭代法是收敛的.（1）为行对角占优阵，即（2)为列对角占优阵，即证明：(1)雅可比迭代矩阵其中（2）为列对角优阵，故为行对角占优阵，由系数矩阵构造的迭代矩阵为行对角占优阵，则有又得到而,得由系数矩阵构造的雅可比迭代矩阵收敛。（如矩阵既是行对角占优阵，也是列对角占优阵)定理4.4若方程组系数矩阵 为对称正定阵，并且也为对称正定，则雅可比迭代收敛。4.2 高斯-塞德尔（GaussSeidel）迭代法高斯-赛德尔迭代计算在雅可比迭代中,用的值代入方程（4.2）中计算出的值，的计算公式是事实上，在计算前，已经得到的值,不妨将已算出的分量直接代入迭代式中，及时使用最新计算出的分量值.因此的计算公式可改为：即用向量计算出的值，用向量计算出的值，用向量计算出的值,这种迭代格式称为高斯塞德尔迭代。对于方程组AX=y ，如果由它构造高斯塞德尔迭代和雅可比迭代都收敛，那么,多数情况下高斯塞德尔迭代比雅可比迭代的收敛效果要好，但是情况并非总是如此。构造方程组的高斯塞德尔迭代格式的步骤与雅可比类似，设将式（4.1）中每个方程的留在方程的左边，其余各项都移到方程的右边；方程两边除以，得到下列同解方程组:记，对方程组对角线以上的取第步迭代的数值，对角线以下的取第步迭代的数值,构造高斯-塞德尔迭代形式：(4.5）例4。2 用高斯塞德尔方法解方程组：解：方程的迭代格式:取初始值有时，时，计算结果如表4。2所示.表4。2 计算结果012340002。52。11.142。50.88 2.0040。98761。621。0042 1.9984 1.00060。12421.0005 2。0002 1.00000.0037高斯塞德尔迭代矩阵设写成等价矩阵表达式：构造迭代形式：有 (4.6)则高斯塞德尔迭代式（4。4）为 （4。7） 称为高斯-塞德尔迭代矩阵。高斯-塞德尔迭代算法高斯-塞德尔迭代的程序实现与雅可比迭代步骤大致相同,对于方程组,在前面的雅可比算法中，假定雅可比迭代矩阵为表示表示,其迭代核心部分是。下面的算法给出由和计算的过程，省略了形成迭代矩阵和对初始化的部分。雅可比迭代的核心部分：WHILEfor(u:=1;u=n,u+） x1u:=x2ufor(i：=1;j=n；j+) s：=gi； for（j：=1；i=n；i+） s:=s+bix2j /注意x2jENDWHILE在高斯赛德尔迭代计算中并不需要形成迭代矩阵，由式（4。5）可知在计算中只要形成矩阵在程序中令向量。它的核心部分是计算迭代式,计算中只需及时将放到的位置上。高斯-塞德尔迭代的核心部分：WHILE for(u:=1;u=n；u+） x1u:=x2u for（i:=1;j=n;j+) s：=gi； for(j：=1；j=n;j+） s：=s+bix2 j /注意x2jx2i：=s ENDWHILE上列算法是在假定迭代收敛的前提下，使用当型(WHILE）结构控制循环。更一般地，可将上列算法中WHILE循环改为FOR循环，通过控制循环次数和观测计算误差终止循环.届将上列算法中WHILE语句改为WHILE 循环次数这时在程序中要增加循环变量的设定和运算.判断高斯塞德尔迭代收敛的方法与判断雅可比迭代收敛类似，一方面从高斯-塞德尔迭代矩阵获取信息,当或的某种范数时，迭代收敛；另一方面，直接根据方程组系数矩阵的特点作出判断。定理4。5 若方程组系数矩阵A为列或行对角优时，则高斯塞德尔迭代收敛。定理4.6 若方程组系数矩阵A为对称正定阵，则高斯塞德尔迭代收敛。例4。3 方程组中,证明当 时GaussSeidel法收敛，而Jacobi迭代法只在时才收敛。解：对法,因为是对称矩阵,因此只要证时正定即可，由顺序主子得而 得于是得到时故对称正定，法收敛。对 法，可根据定理4.2，由于迭代矩阵 即 是 法收敛的充要条件，故法只在时才收敛。当时，法收敛，而，法不收敛，此时不是正定的.4.3 逐次超松弛（SOR）迭代法逐次超松弛迭代计算方程组的雅可比迭代形式记其中是下三角矩阵，是上三角矩阵。得高斯-塞德尔的迭代形式:记，有这样可视为的修正量，而恰是由加修正量而得，如果将改为)加上修正量乘一个因子，迭代格改为：整理得 (4。8)这里为修正因子，称为松弛因子,而式（4.8）称为松弛迭代。迭代的分量形式为(4。9)式(4.9)称为松弛迭代的计算格式。例4.4 给定方程组用SOR法求解，取解：用SOR迭代公式可得取初始值：如果用高斯赛德尔迭代法迭代72次得：用SOR迭代法，只须迭代25次即得: 逐次超松弛迭代算法下列算法假定迭代矩阵收敛，否则要在WHILE循环中增加判断条件.1定义和输入系数矩阵与常数项向量的元素，输入松弛因子的值。2FOR i:=1,2，n / 假定，形成常数项向量 FOR 34。 WHILE ; ENDWHILE5输出.松弛迭代矩阵将式（4.9）中含有与的项分别放在方程两边： 用 代入上式，有 则松弛因子为的迭代矩阵为 其中 定理4。7 逐次超松弛迭代收敛的必要条件。定理4.8 若为正定矩阵，则当时，逐次超松弛迭代恒收敛。以上定理给出了逐次超松弛迭代因子的范围.对于每个给定的方程组，确定究竟取多少为最佳，这是比较困难的问题，对某些特定的方程组，我们可以得到一些理论结果.通常，把的迭代称为亚松弛迭代，把的迭代称为高斯塞德尔迭代，而把的迭代称为松弛迭代。4.4逆矩阵计算在线性代数中逆矩阵是按其伴随矩阵定义的，若则方阵可逆，且，其中为的伴随矩阵。要计算个阶的列式才能得到一个伴随矩阵，在数值计算中因其计算工作量大而不被采用。通常对做行的初等的效换，在将化成的过程中得到.在数值计算中，这仍然是一种行之有效的方法。由逆矩阵的定义令，有化为个方程组j是第个分量为1，其余分量为0的维向量。或记为:。用直接法或迭代法算出也就完成了逆矩阵计算。如果依次对用高斯若尔当消元法，组合起来看有（当然也能组合起来做）：这正是在线性代数中用初等变换计算逆矩阵的方法.由此可见,计算一个阶逆矩阵的工作量相当于解个线性方程组。在数值计算中常常将计算矩阵逆的问题转化为解线性方程组的问题。例如，已知方阵和向量有迭代关系式，在计算中不是先算出，再作与的乘积得到；而将作为线性方程组系数矩阵,求解方程组作为常驻数项解出.4.5程序示例程序4。1用雅可比迭代求解方程组：算法描述输入系数矩阵及常数项向量c。按雅可比迭代公式：求解。程序源码/ Purpose:雅可比迭代求解线性方程组 /includeincludestdio.hdefine MAXN 20 /方程的最大维数define MAXREPT 100#define epsilon 0.00001 /求解精度int main（ ） int n; int I, j, k; double err； static double a MAXN MAX-N,bMAXN MAX-N，cMAX-N，g MAX-N； static double xMAXN,nxMAXN； printf（nInput n value(dim of Ax=c）：)； /输入方程的维数 scanf(%d,&n)； if（nMAX-N) printf（The input n is larger then MAX-N,please redefine the MAX-N.n）； if（n=0) printf(please input a number between 1 and %d.n), MAXN；return 1;/输入Ax=c的系数矩阵A.PRINTF（Now input the matrix a（I，j，)=0，%d：n，n-1)；for(i=0; in; i+) /形成迭代矩阵b for(j=0; jn； j+)bij=aij/aii；gi=ci/aij; /为了简化程序，假设aii！0,/否则要附加对aii的判断及其相应的处理for（i = 0； i MAXREPT; i+） for(j=0; jn; j+ ） nxj=gj； for(j=0, jn; j+ ）for（k=0, kn; k+ )if (j=k）continue;nxj+ =bjkxk； /迭代 err = 0 ；for（j=0, jn; j+ ）if (errfabs（nxj-xj）err=fabs（nxj-xj）； /误差for（j=0， jn； j+ )xj=nxj；if(errepsilon)printf（“Solvex_I=n”）; /输出for(i=0， in； i+ )printf(“%fn”，xi）;return 0； printf（“After% d repeat, no result n， MAXREPT ); /输出return 1; 计算实例用雅可比方法解下列方程组： 程序输入输出Inprt n value(dim of Ax=c）：3Now input the matrix a(j, j）, i， j =0，2:64 13 1 2 90 1 1 1 40Now input the matrix c(i)， i= 0，, 2：14 -5 20Solvex_i=0。2295470。0661300.492608程序4。2用逐次超松弛迭代（ 作参数)求解方程组： 算法描述输入矩阵及列向量c.按因子为的超松弛迭代公式：求解.程序源码/Purpose：超松弛迭代求解线性方程组 /# include# includeMAX_N） printf(“The input nis larger then MAX_N, please redefine the MAX_N。n”); return 1;if （n=0) printf(“Please input a number between 1 and % d。 n”， MAX_N)； return 1； / /输入的A矩阵printf（“Now input the matrix a（i, j）， i, j =0，% d： n， n-1);for（i=0； in; i+)for(j=0； jn; j+） scanf(“%1f” ， &sij)；printf（“Now input the matrix c（i）， i =0， d: n, n1）;for(i=0; in； i+) scanf （“1f” ， &ci)；printf(“Now input the w value:）;scanf(“%1f ， w）;if (w1| w=2）printf(“w must between 1and 2。 n”);return 1;for（i=0; in； i+） /改造x k+1=bx_k+g迭代矩阵for（j=0; jn; j+） bij=-aij/aii；gi=ci/aii；for(i=0； iMAXREPT； i+）for(j=0； jn； j+）nxj=gj;for(j=0; jn; j+)for(k=0； kj; k+）err=0;for(j=0; jn ； j+）if（errfabs(nxj-xj)）err=fabs(nxjxj); /误差 for(j=0; jn ; j+)xj=nxj；if（errepsilon)printf（“Solve x_i = n”) ； /输出 for(i=0； in ; i+)printf(“fn”,xi）;return 0；printf（“After % d repeat, no result n， MAXREPT）； /输出return 1 ；计算实例解下列方程组：程序输入输出Input n value(dim of Ax=c）: 3Now input the matrix a(i， j)， i， j=0，2：64 3 -1 2 -90 1 1 1 40Now input the matrix c（i)， i=0,，2：14 -5 20Now input the w value： 1 Solve x_i=0.2295470。0661300.492608