首页 > 算法 > 正文

阅读排行

宋叶志书 - 选主元消去法之我见
2014-04-02 12:21:59   来源:岸边的鱼 @ Fcode研讨团队   评论:0 点击:

针对宋叶志老师编著的一书中,第一章第三节中列主元消去法不能正确选取出列主元的问题,进行举例说明并给出改变方案。

写在前面的话:感谢fcoder群友一直以来对我的无私帮助,才有我今天对fortran语言的热爱,感谢群主大哥的建议才有本文的诞生。

最近在使用宋叶志[1]老师编著的《FORTRAN95/2003科学计算与工程》一书中的第一章,第三节选主元消去法的程序时,偶然发现程序选取出的主元老是有问题,遂使用debug一步步调试终于发现,原来是宋老师的程序有点小小的问题,现斗胆写出来,请广大群友批评指正。

在求解线性方程组时,通常的直接解法包括,高斯消去法,高斯选主元消去法,LU分解法等多种方法。

个人认为以高斯消去法的思想最朴素,也最容易被初学者所接受。该方法通过对增广矩阵实施消元变换,而在变换的过程中与原方程始终保持等价,直到矩阵变为上三角矩阵,再采用解三角方程组的方法进行回代,便得到原方程组的解。

然而,这种朴素的思想当中蕴涵着一个致命的缺憾,即如果方程的增广矩阵在消元的过程中出现0主元或者主元非常小的话,高斯消去法将失败或者数值不稳定。为了克服这一缺憾,发展出了高斯选主元消去法(包括列主元消去法,行主元消去法,全主元消去法三种)。下面给出列主元消去法的算法:
(1)设置增广矩阵A_{p}=\begin{bmatrix}A,b\end{bmatrix}
(2)对K=1,N-1,做(3)---(6)处理。
(3)设置一个元素a_{max}=\left | A_{p}(k,k) \right |,及标号ID_{max}=k
(4)查找当列的最大元素,然后用标号ID_{max}记录下这个元素所在的行数。
(5)交换第k行与第ID_{max}行的所有数据,注意其他元素保持不变。
(6)完成(5)时已经完成了主元的选取,下面可以按照消去法进行计算。
(7)完成以上步骤之后,已经形成上三角矩阵。
(8)对上三角矩阵进行回代,即可得到方程的解。
   下面我们再来重点看下计算过程中的第(4)步:查找当列的最大元素,这一说法到底应该怎么理解?是查找出来的元素最大,还是查找出来的元素的绝对值最大?笔者结合第(3)步中我们使用的初始比较值是\left | A_{p}(k,k) \right |,认为应该理解为查找出来的元素的绝对值最大,宋叶志老师书中给出的选取主元的程序是:

do k=1,N-1
  elmax=dabs(Ab(k,k))
  id_max=k
  !这段为查找主元素	
  !这段程序的主要目的不是为了赋值最大元素给 elmax 而是为了找出最大元素对应的标号
	do i=k+1,n
    if (dabs(Ab(i,k))>elmax) then
      elmax=Ab(i,k)
      id_max=i     
    end if          
  end do  
!至此,已经完成查找最大元素,查找完成以后与第k行交换

现在我们按照宋老师的程序来看下面的例子:假定第一列元素为\left [ 1,2,-5,1 \right ]^{-1},(此时n=4)我们选取的a_{max}=\left | A_{p}(k,k) \right |=\left | A_{p}(1,1) \right |=1,此时k=1,elmax=1,当i=2时开始执行if比较:
第1次i=2,dabs(Ab(i,k))>elmax即2>1成立,此时:elmax=2,id_max=2。
第2次i=3,dabs(Ab(i,k))>elmax即5>2成立,此时:elmax=-5,id_max=3。
第3次i=4,dabs(Ab(i,k))>elmax即1>-5成立,此时:elmax=1,id_max=4。
第4次i=5,结束if判断。

从上面的执行情况来看,按照宋老师的程序,我们选出来的主元将会是1,而1既不是此列中的最大元素,也不是绝对值最大的元素,即选取主元失败。

分析原因不难发现失败的原因在于列中含有-5这个元素,它的存在使得程序选取出的主元既不是最大的,也非绝对值最大的。当然程序的修订也很简单,我们只需要将elmax=Ab(i,k),修改为elmax=dabs(Ab(i,k))即可。不难发现,修改之后程序即可正确选择出该列中的绝对值最大的元素-5。

参考文献

[1] 宋叶志, 茅永兴, 赵秀杰. Fortran95/2003科学计算与工程[M]. 北京:清华大学出版社, 2011. 12-16

相关热词搜索:高斯列主元消去法

上一篇:Fortran求大数阶乘
下一篇:高斯勒让德求积分Fortran程序

分享到: 收藏