最近在使用宋叶志[1]老师编著的《FORTRAN95/2003科学计算与工程》
在求解线性方程组时,通常的直接解法包括,高斯消去法,高斯选主元消去法,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