第二章 线性方程组的直接方法(1)
2.1 顺序高斯消去法
很多工程实际问题最终都可以归结为求解线性方程组。如:
-
用最小二乘法求实验数据的曲线拟合问题
-
工程中的三次样条函数的插值问题
-
经济运行中的投入产出问题
-
大地测量、机械与建筑结构的设计计算问题
这些问题都可以归结为求解线性方程组或者非线性方程组的数学问题。
本节讨论n元线性方程组的直接解法
⎩⎪⎪⎪⎨⎪⎪⎪⎧a11x1+a12x2+...+a1nxn=b1a21x1+a22x2+...+a2nxn=b2...............................an1x1+an2x2+...+annxn=bn
矩阵形式为Ax=b,为
A=⎣⎢⎢⎡a11a21...an1a12a22...an2a13.........a1na2n...ann⎦⎥⎥⎤,x=⎣⎢⎢⎡x1x2...xn⎦⎥⎥⎤,b=⎣⎢⎢⎡b1b2...bn⎦⎥⎥⎤
根据线性代数知识我们知道,这个矩阵有唯一解的条件是det(A)=0
求解线性方程组的直接方法:如果不考虑计算过程当中的舍入误差,经过有限次算数运算就可以得到方程组的精确解。
但是我们知道,我们的研究方法都要拿到计算机上实现,舍入误差是不可避免地。上讲提到的Cramer法则就是一种直接方法,但是它并不实用。如果用Cramer法则去计算线性方程组的话,计算量非常大。
接下来介绍几种实用的直接方法。
Gauss消去法
它是一种规则化的加减消元法,它的基本思想是通过逐次消元计算,把一个一般的线性方程组的求解问题转化为等效的上三角方程组的求解问题。
简单示例
⎩⎪⎨⎪⎧2x1+4x2−2x3=2x1−3x2−3x3=−14x1+2x2+2x3=3
消去后面两式x1得
⎩⎪⎨⎪⎧2x1+4x2−2x3=2−5x2−2x3=−2−6x2+6x3=−1
消去第三个方程中的x2,得
⎩⎪⎨⎪⎧2x1+4x2−2x3=2−5x2−2x3=−2542x3=57
这个时候方程已经变成一个上三角形式的方程组,所以我们通过回代求解出x1,x2,x3,得
x3=61,x2=31,x1=21
这就是顺序高斯消元法的基本思想,先逐次消元,然后回代求解。
下面用矩阵形式来将刚刚的消元过程表示,需要注意的是,上面的消元过程是对增广矩阵来进行的,而不是系数矩阵。
(A,b)=⎣⎡2144−32−2−322−13⎦⎤
消去后面两行中的第一列,r2−21r1,r3−2r1得
⎣⎡2004−5−6−2−262−2−1⎦⎤
消去最后一行的第二列,r3−56r2得
⎣⎡2004−50−2−25422−257⎦⎤
消元过程实际上就是对增广矩阵不断地实施初等行变换。
一般方法
接下来来介绍求解线性方程组(1)的顺序Gauss消去法,为了清楚起见,我们记A(1)=A,b(1)=b,aij1=aij,bi(1)=bi
则,线性方程组(1)的增广矩阵为:
(A(1),b(1))=⎣⎢⎢⎢⎢⎢⎡a11(1)a21(1)a31(1)...an1(1)a12(1)a22(1)a32(1)...an2(1)a13(1)a23(1)a33(1)...an3(1)...............a1n(1)a2n(1)a3n(1)...ann(1)b1(1)b2(1)b3(1)...bn(1)⎦⎥⎥⎥⎥⎥⎤
第一步:设a11(1)=0,依次用−li1=−a11(1)ai1(1),(i=2,3,...,n)乘矩阵的第1行加到第i行,得到矩阵:
(A(2),b(2))=⎣⎢⎢⎢⎢⎢⎡a11(1)00...0a12(1)a22(2)a32(2)...an2(2)a13(1)a23(2)a33(2)...an3(2)...............a1n(1)a2n(2)a3n(2)...ann(2)b1(1)b2(2)b3(2)bn(2)⎦⎥⎥⎥⎥⎥⎤
其中:aij(2)=aij(1)−li1a1j(1),i,j=2,3,...,nbi(2)=bi(1)−li1b1(1),i=2,3,...,n
第二步:设a22(2)=0,依次用−li2=−a22(2)ai2(2),(i=3,4,...,n)乘矩阵的第2行加到第i行,得到矩阵:
(A(3),b(3))=⎣⎢⎢⎢⎢⎢⎡a11(1)00...0a12(1)a22(2)0...0a13(1)a23(2)a33(3)...an3(3)...............a1n(1)a2n(2)a3n(3)...ann(3)b1(1)b2(2)b3(3)...bn(3)⎦⎥⎥⎥⎥⎥⎤
其中:aij(3)=aij(2)−li2a2j(2),i,j=3,4,...,nbi(3)=bi(2)−li2b2(2),i=3,4,...,n
如此继续消元下去,第n−1步结束后得到矩阵:
(A(n),b(n))=⎣⎢⎢⎢⎢⎢⎡a11(1)00...0a12(1)a22(2)0...0a13(1)a23(2)a33(3)...0...............a1n(1)a2n(2)a3n(3)...ann(n)b1(1)b2(2)b3(3)...bn(n)⎦⎥⎥⎥⎥⎥⎤
观察这个矩阵,第一行矩阵从头到尾均没有发生变化,第二行元素只发生一次变化,第三行发生了两次变化,这里我们发现,系数矩阵A已经变成了上三角矩阵U,说明消元过程已经结束了。
对应的方程变成:
⎩⎪⎪⎪⎨⎪⎪⎪⎧a11(1)x1+a12(1)x2+...+a1n(1)xn=b1(1)a22(2)x2+...+a2n(2)xn=b2(2)...............................ann(n)xn=bn(n)
对方程组进行回代,就可求出方程组的解。
整个消元和回代的过程加在一起就是我们所说的顺序Gauss消去法,那么大家看到了,在整个过程中,方法是精确的,也就是没有截断误差,也就是说顺序Gauss消去法确实是一种直接方法。那么它对比于Cramer法则,它的计算量是如何呢?
算法效率
顺序Gauss消去法求解n元线性方程组的乘除运算量是:
n2−1+(n−1)2−1+...+22−1+1+2+...+n=∑k=1n(k2−1)+∑k=1nk=6n(n+1)(2n+1)−n+2n(n+1)=31(n3+3n2−n)
例如,假设n=20,则乘除次数为3060次,Cramer法则为9.7×1020
因此,顺序Gauss消去法确实大大减少了运算量,因此它是一种比较实用的直接方法,简称为Gauss消去法。
局限性
(1)主元为0造成的局限性
需要注意的是,在消元之前,我们假设了akkk不能为0,如果为0则无法进行消元,我们将akkk(k=1,2,...,n)称为主元素。这里的k需要真正算到这一步的时候才能判断出k是否为0,如果我们在某一步的时候才判断出k=0了,那么前面的步骤就全部白算了。因此我们希望有一种判断方法,能在开始消元之前就判断出整个Gauss消去法中,主元素能否为0。结论为:
主元素都不为0的充要条件是:矩阵A的各阶顺序主子式都不为0
矩阵A为系数矩阵,因此我们可以根据判断矩阵A的各阶顺序主子式来判断,因此这也是Gauss消去法的局限所在。
(2)消去过程中产生“大数吃小数”问题
顺序高斯消去法除了当主元素为0时可能造成的无法消元以外还有其他局限性吗?我们来看一个例子
例1:解线性方程组(用十进制四位浮点计算)
{0.000100x1+1.00x2=1.001.00x1+1.00x2=2.00
(用Cramer法则可得精确解x1∗=1.00010,x2∗=0.99990)
解:用顺序Gauss消去法,消元得
{0.000100x1+1.00x2=1.00−10000x2=−10000
这就产生了“大数吃小数”的问题,这样求解得x2=1.00,x1=0.00。我们发现这个误差非常非常大。由于顺序Gauss消去法没有截断误差,因此这是由舍入误差造成的。那么这个“大数”是从哪来的呢?我们发现在第一次进行消元的时候,求两行比例系数的过程中,我们用了绝对值很小的数做了除数,以至于误差被放大了。
我们将方程组交换两行,改写成
{1.00x1+1.00x2=2.000.000100x1+1.00x2=1.00
用顺序Gauss消去法,消元得:
{1.00x1+1.00x2=2.001.00x2=1.00
回代得解x2=1.00,x1=1.00,这个近似解是可以接受的,误差也在可接受的范围。
通过这个例子大家看到了,在顺序Gauss消去法当中,尽管主元素不为0,可以保证Gauss消去法从头进行到尾,但是如果在消元的过程中,主元素很小的话,同样会产生最终解的舍入误差非常大以至于不能用的情况。这就是顺序高斯消去法第二个局限性,同时我们也找到了改变这个局限性的方法,就是互换两个方程。