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

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

———————————————————————————————— 作者：
———————————————————————————————— 日期：

30

个人收集整理 勿做商业用途
第4章  解线性方程组的迭代法
　　用迭代法求解线性方程组与第4章非线性方程求根的方法相似，对方程组进行等价变换，构造同解方程组(对可构造各种等价方程组，如分解，可逆，则由得到），以此构造迭代关系式
　　　　　　　　　　    　          　 （4。1）
　　任取初始向量，代入迭代式中，经计算得到迭代序列。
　　若迭代序列收敛，设的极限为，对迭代式两边取极限

　　即是方程组的解,此时称迭代法收敛，否则称迭代法发散。我们将看到，不同于非线性方程的迭代方法，解线性方程组的迭代收敛与否完全决定于迭代矩阵的性质，与迭代初始值的选取无关。迭代法的优点是占有存储空间少,程序实现简单,尤其适用于大型稀疏矩阵；不尽人意之处是要面对判断迭代是否收敛和收敛速度的问题.
　　可以证明迭代矩阵的与谱半径是迭代收敛的充分必要条件，其中是矩阵的特征根。事实上，若为方程组的解，则有

　　再由迭代式可得到

　　由线性代数定理，的充分必要条件.
　　因此对迭代法（4。1）的收敛性有以下两个定理成立。
　　定理4.1 迭代法 收敛的充要条件是。
　　定理4.2 迭代法收敛的充要条件是迭代矩阵的谱半径
　　因此,称谱半径小于1的矩阵为收敛矩阵。计算矩阵的谱半径，需要求解矩阵的特征值才能得到，通常这是较为繁重的工作.但是可以通过计算矩阵的范数等方法简化判断收敛的工作。前面已经提到过，若|｜A｜｜p矩阵的范数，则总有。因此，若，则必为收敛矩阵。计算矩阵的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．定义和输入系数矩阵与常数项向量的元素。
　　2．FOR i：=1,2,…，n
　　　　｛   //假定，形成常数项向量
　　            FOR j:=1，2,…，n

　　｝              //形成迭代矩阵元素
　　3．
　　                  // 赋初始值，x1和x2分别表示和
　　4．WHILE   
　　     x1：=x2
　　     x2：=B＊x1+g   //  FOR  u：=1，2，…,n
　　                  //     s:= g［u］；
　　                  //  FOR v:=1，2,…，n   s：=s+b[u]［v]＊x1［v］；
　　                  //     x2［u］:=s；
　　  ENDWHILE
　　5．输出方程组的解
　　例4.1 用雅可比方法解下列方程组:

　　解：方程的迭代格式：

　　或

　　雅可比迭代收敛.
　　取初始值,计算结果由表4。1所示。
表4.1 计算结果

0
　　1
　　1
　　1

1
　　—1。5
　　1。6
　　0.9
　　0.25
2
　　—1.25
　　2。08
　　1.09
　　0.48
3
　　-0。915
　　2.068
　　1。017
　　0。335
4
　　—0.9575
　　1.9864
　　0.9847
　　0。0816
5
　　—1。01445
　　1。98844
　　0.99711
　　0.05695
6
　　-1。00722
　　2。00231
　　1.0026
　　0.01387
7
　　—0.997543
　　2.00197
　　1。00049
　　0。009687
　　方程组的准确解是
　　4。1.2 雅可比迭代收敛条件
　　对于方程组，构造雅可比迭代格式其中，。当迭代矩阵的谱半径时，迭代收敛，这是收敛的充分必要条件。迭代矩阵的某范数时，迭代收敛。要注意的是范数小于1只是判断迭代矩阵收敛的充分条件，当迭代矩阵的一种范数||B|｜〉1，并不能确定迭代矩阵是收敛还是发散。例如，，则，但它的特征值是0.9和0.8。是收敛矩阵。
　　当方程组的系数矩阵具有某些性质时，可直接判定由它生成的雅可比迭代矩阵是收敛的.
　　定理4。3 若方程组的系数矩阵，满足下列条件之一，则其雅可比迭代法是收敛的.
　　（1）为行对角占优阵，即
　　（2)为列对角占优阵，即
　　证明：(1)雅可比迭代矩阵其中

　　（2）为列对角优阵，故为行对角占优阵，由系数矩阵构造的迭代矩阵为行对角占优阵，则有

　　又
　　得到
　　而,
　　得
　　由系数矩阵构造的雅可比迭代矩阵收敛。
　　（如矩阵既是行对角占优阵，也是列对角占优阵)
　　定理4.4　若方程组系数矩阵 为对称正定阵，并且也为对称正定，则雅可比迭代收敛。

4.2  高斯-塞德尔（Gauss—Seidel）迭代法
　　高斯-赛德尔迭代计算
　　在雅可比迭代中,用的值代入方程（4.2）中计算出的值，的计算公式是

　　事实上，在计算前，已经得到的值,不妨将已算出的分量直接代入迭代式中，及时使用最新计算出的分量值.因此的计算公式可改为：

　　即用向量计算出的值，用向量计算出的值，用向量计算出的值,这种迭代格式称为高斯—塞德尔迭代。
　　对于方程组AX=y ，如果由它构造高斯—塞德尔迭代和雅可比迭代都收敛，那么,多数情况下高斯—塞德尔迭代比雅可比迭代的收敛效果要好，但是情况并非总是如此。
　　构造方程组的高斯—塞德尔迭代格式的步骤与雅可比类似，设将式（4.1）中每个方程的留在方程的左边，其余各项都移到方程的右边；方程两边除以，得到下列同解方程组:

　　记，对方程组对角线以上的取第步迭代的数值，对角线以下的取第步迭代的数值,构造高斯-塞德尔迭代形式：

(4.5）
　　例4。2 用高斯—塞德尔方法解方程组：

　　解：方程的迭代格式:

　　取初始值有
　　时，

　　时，

　　计算结果如表4。2所示.
表4。2 计算结果

0
1
2
3
4
　　　0　　　　　0　　　　　0
　　　－2。5　　　2。1　　　　1.14　　　　　2。5
　　　－0.88　　 2.004　　　0。9876　　　　1。62
　　　－1。0042　 1.9984　　 1.0006　　　　0。1242
　　　－1.0005　 2。0002　　 1.0000　　　　0.0037
　　高斯—塞德尔迭代矩阵
　　设

　　写成等价矩阵表达式：

　　构造迭代形式：

　　有      　　　　(4.6)
　　则高斯—塞德尔迭代式（4。4）为
　　　　　　　　          　　　　　   （4。7）

　　称为高斯-塞德尔迭代矩阵。
　　高斯-塞德尔迭代算法
　　高斯-塞德尔迭代的程序实现与雅可比迭代步骤大致相同,对于方程组,在前面的雅可比算法中，假定雅可比迭代矩阵为表示表示,其迭代核心部分是。下面的算法给出由和计算的过程，省略了形成迭代矩阵和对初始化的部分。
　　雅可比迭代的核心部分：
　　WHILE
　　for(u:=1;u<=n,u++）  x1［u]:=x2［u]
　　for(i：=1;j〈=n；j++)
　　｛   s：=g[i］；
　　       for（j：=1；i<=n；i++）
　　              ｛ s:=s+b［i]x2[j]} //注意x2[j］

　　｝
　　ENDWHILE
　　在高斯—赛德尔迭代计算中并不需要形成迭代矩阵，由式（4。5）可知在计算中只要形成矩阵在程序中令向量
　　。它的核心部分是计算迭代式,计算中只需及时将放到的位置上。
　　高斯-塞德尔迭代的核心部分：
　　WHILE    
　　         for(u:=1;u〈=n；u++）  x1［u］:=x2［u]
　　           for（i:=1;j〈=n;j++)
　　                 ｛ s：=g[i]；
　　              for(j：=1；j〈=n;j++）
　　                ｛ s：=s+b［i］x2［ j ］}     //注意x2［j］
　　                  x2[i]：=s }
　　ENDWHILE
　　上列算法是在假定迭代收敛的前提下，使用当型(WHILE）结构控制循环。更一般地，可将上列算法中WHILE循环改为FOR循环，通过控制循环次数和观测计算误差终止循环.届将上列算法中WHILE语句改为
WHILE    循环次数
　　这时在程序中要增加循环变量的设定和运算.
　　判断高斯塞德尔迭代收敛的方法与判断雅可比迭代收敛类似，一方面从高斯-塞德尔迭代矩阵获取信息,当或的某种范数时，迭代收敛；另一方面，直接根据方程组系数矩阵的特点作出判断。
　　定理4。5 若方程组系数矩阵A为列或行对角优时，则高斯塞德尔迭代收敛。
　　定理4.6 若方程组系数矩阵A为对称正定阵，则高斯塞德尔迭代收敛。
　　例4。3 方程组中,,
　　证明当 时Gauss—Seidel法收敛，而Jacobi迭代法只在时才收敛。
　　解：对法,因为是对称矩阵,因此只要证时正定即可，由顺序主子
　　得
　　而 得于是得到时故对称正定，法收敛。
　　对 法，可根据定理4.2，由于迭代矩阵

　　即
　　是 法收敛的充要条件，故法只在时才收敛。
　　当时，法收敛，而，法不收敛，此时不是正定的.

4.3  逐次超松弛（SOR）迭代法
　　逐次超松弛迭代计算
　　方程组的雅可比迭代形式记其中是下三角矩阵，是上三角矩阵。
　　得高斯-塞德尔的迭代形式:

　　记，有

　　这样可视为的修正量，而恰是由加修正量而得，如果将改为)加上修正量乘一个因子，迭代格改为：

　　整理得
　　　　　     (4。8)
　　这里为修正因子，称为松弛因子,而式（4.8）称为松弛迭代。迭代的分量形式为
(4。9)
　　式(4.9)称为松弛迭代的计算格式。
　　例4.4 给定方程组

　　用SOR法求解，取
　　解：用SOR迭代公式可得

　　取初始值：
　　如果用高斯—赛德尔迭代法迭代72次得：
　　用SOR迭代法，只须迭代25次即得:
　　逐次超松弛迭代算法
　　下列算法假定迭代矩阵收敛，否则要在WHILE循环中增加判断条件.
　　1．定义和输入系数矩阵与常数项向量的元素，输入松弛因子的值。
　　2．FOR i:=1,2，…，n
　　        // 假定 ，形成常数项向量
　　              FOR

　　       ｝
　　3．
　　4。  WHILE 

　　                ;

　　                      ｝
　　　ENDWHILE
　　5．输出.
　　松弛迭代矩阵
　　将式（4.9）中含有与的项分别放在方程两边：

　　用   代入上式，有

　　则松弛因子为的迭代矩阵为

　　其中        

　　定理4。7 逐次超松弛迭代收敛的必要条件。
　　定理4.8 若为正定矩阵，则当时，逐次超松弛迭代恒收敛。
　　以上定理给出了逐次超松弛迭代因子的范围.对于每个给定的方程组，确定究竟取多少为最佳，这是比较困难的问题，对某些特定的方程组，我们可以得到一些理论结果.
　　通常，把的迭代称为亚松弛迭代，把的迭代称为高斯—塞德尔迭代，而把的迭代称为松弛迭代。
4.4  逆矩阵计算
　　在线性代数中逆矩阵是按其伴随矩阵定义的，若则方阵可逆，且，其中为的伴随矩阵。要计算个阶的列式才能得到一个伴随矩阵，在数值计算中因其计算工作量大而不被采用。通常对做行的初等的效换，在将化成的过程中得到.在数值计算中，这仍然是一种行之有效的方法。
　　由逆矩阵的定义 令，有

　　化为个方程组
j
　　是第个分量为1，其余分量为0的维向量。或记为:。
　　用直接法或迭代法算出也就完成了逆矩阵计算。
　　如果依次对用高斯若尔当消元法，组合起来看有（当然也能组合起来做）：

　　这正是在线性代数中用初等变换计算逆矩阵的方法.
　　由此可见,计算一个阶逆矩阵的工作量相当于解个线性方程组。在数值计算中常常将计算矩阵逆的问题转化为解线性方程组的问题。
　　例如，已知方阵和向量有迭代关系式，在计算中不是先算出，再作与的乘积得到；而将作为线性方程组系数矩阵,求解方程组作为常驻数项解出.

4.5　程序示例
　　程序4。1　用雅可比迭代求解方程组：

　　算法描述
　　输入系数矩阵及常数项向量c。
　　按雅可比迭代公式：

　　求解。
　　程序源码
　　/////////////////////////////////////////
　　// 　Purpose:雅可比迭代求解线性方程组  //
　　/////////////////////////////////////////
　　＃include<stdio。h>
　　＃include〈stdio.h>
　　＃define MAX—N 20         //方程的最大维数
　　＃define MAXREPT 100
　　#define epsilon 0.00001      //求解精度
　　int main（ ）
　　｛  int n;
　　   int I, j, k;
　　   double err；
　　   static double a ［MAX—N］ ［MAX-N],b［MAX—N] ［MAX-N］，c[MAX-N]，g [MAX-N]；
　　   static double x［MAX—N］,nx［MAX—N]；
　　   printf（〃\nInput n value(dim of Ax=c）：〃)；  //输入方程的维数
　　   scanf(〃%d〃,&n)；
　　   if（n>MAX-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)〃, MAX—N｝；return 1;｝
　　//输入Ax=c的系数矩阵A.
　　PRINTF（〃Now input the matrix a（I，j，)=0，…，%d：\n〃，n-1)；
　　for(i=0; i<n; i++)   //形成迭代矩阵b
　　  for(j=0; j<n； j++)    
　　{
　　b[i]［j]=－a[i］[j］/a［i]［i］；
　　g[i］=c[i]/a［i][j］;        //为了简化程序，假设a［i］［i]！≠0,
　　//否则要附加对a[i］[i]的判断及其相应的处理
　　{
　　for（i = 0； i< MAXREPT; i++）
　　{ for(j=0; j〈n; j++ ）
　　  nx[j]=g[j]；
　　  for(j=0, j〈n; j++ ）
　　｛
　　for（k=0, k〈n; k++ )
　　{
　　if (j==k）continue;
　　nx[j]+ =b[j］［k］＊x［k］；    //迭代
　　}
　　  ｝
　　err = 0 ；
　　for（j=0, j<n; j++ ）
　　if (err〈fabs（nx[j］-x［j]））err=fabs（nx[j］-x[j］）；   //误差
　　for（j=0， j<n； j++ )
　　x［j]=nx[j］；
　　if(err<epsilon)
　　｛
　　printf（“Solve…x_I=\n”）;     //输出
　　for(i=0， i〈n； i++ )printf(“%f\n”，x［i]）;
　　return 0；
　　  {
　　｝
　　printf（“After% d repeat, no result…\ n"， MAXREPT );    /输出
　　return 1; {
　　计算实例
　　用雅可比方法解下列方程组：

　　程序输入输出
　　Inprt n value(dim of Ax=c）：3
　　Now input the matrix a(j, j）, i， j =0，…2:
　　64  —13  —1  2  –90  1  1  1  40
　　Now input the matrix c(i)， i= 0，…, 2：
　　14  -5  20 
　　Solve…x_i=
　　0。229547
　　0。066130
　　0.492608
　　程序4。2　用逐次超松弛迭代（ 作参数)求解方程组：

　　算法描述
　　输入矩阵及列向量c.
　　按因子为的超松弛迭代公式：

　　 求解.
　　程序源码
　　/////////////////////////////////////////
　　//    Purpose：超松弛迭代求解线性方程组 //
　　/////////////////////////////////////////
　　# include<stdio.h>
　　# include<math。h〉
　　# define MAX_N 20      / /方程的最大维数
　　＃ define MAXREPT 100
　　＃ define epsilon 0。00001      / /求解精度
　　int main（ ）
　　{ int n;
　　int i, j， k ;
　　double err, w;
　　static double a［MAX_N］[MAX_N］， b［MAX_N］［MAX_N]，c［MAX_N]，c［MAX_N]，g［MAX_N] ;
　　static double a[MAX_N］ , nx［MAX_N］;
　　PRINTF（“\nlnput n value（dim of Ax=c ）：”)； //输入方程的维数
　　Scanf（“% d”, &n)；
　　If （n>MAX_N）
　　{  printf(“The input n—is 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； i〈n; i++)
　　for(j=0； j〈n; j++）
　　  scanf(“%1f” ， &s[i］[j］)；
　　printf（“Now input the matrix c（i）， i =0，…％ d: \n", n—1）;
　　for(i=0; i〈n； i++) scanf （“％1f” ， &c[i］)；
　　printf(“Now input the w value:"）;
　　scanf(“%1f" ， ＆w）;
　　if (w〈1|| w>=2）
　　}
　　printf(“w must between 1and 2。 \ n”);
　　return 1;
　　｛
　　for（i=0; i〈n； i++）   //改造x ｛k+1｝=bx_｛k｝+g迭代矩阵
　　for（j=0; j〈n; j++）
　　  {
　　  b［i]［j]=-a[i］[j]/a［i][i］；
　　g[i］=c［i]/a［i］［i]；
　　｛
　　for(i=0； i<MAXREPT； i++）
　　｛
　　for(j=0； j〈n； j++）
　　nx［j］=g［j];
　　for(j=0; j〈n; j++)
　　{
　　for(k=0； k〈j; k++）

　　｛
　　err=0;
　　for(j=0; j<n ； j++）
　　if（err<fabs(nx[j］-x［j］)）err=fabs(nx[j]—x［j]);   //误差
　　  for(j=0; j〈n ; j++)
　　x[j]=nx［j]；
　　if（err<epsilon)
　　{
　　printf（“Solve … x_i = \ n”) ；    //输出
　　   for(i=0； i〈n ; i++)printf(“％f\n”,x［i］）;
　　return 0；
　　}
　　 {
　　 printf（“After % d repeat, no result … \ n"， MAXREPT）；    //输出
　　 return 1 ；
　　计算实例
　　解下列方程组：

　　程序输入输出
　　Input n value(dim of Ax=c）: 3
　　Now input the matrix a(i， j)， i， j=0，…，2：
　　64  —3  -1  2  -90  1  1  1  40
　　Now input the matrix c（i)， i=0,…，2：
　　14  -5  20
　　Now input the w value： 1
　　    Solve… x_i=
　　0.229547
　　0。066130
　　0.492608