消元法

线性方程组的直接法

张瑞
中国科学技术大学数学科学学院

消元法

线性方程组可以写成一般矩阵形式形式

\[\begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \cdots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{pmatrix} = \begin{pmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{pmatrix} \]

简记为

\[Ax=b \]

根据线性代数的知识,当且仅当$\det A\neq 0$时,方程组存在唯一的解。

Gramer法则给出了解的表达式

\[x_i=\frac{\det A_i}{\det A}, i=1,2,\cdots,n \]

其中$A_i$是把$A$的第$i$列换成$b$后得到的方阵。

. 但这个公式不能用于计算机求解,因为计算量太大。

直接法:

  • 经过若干次等价变换后再求解,因而该类方法没有截断误差。
  • 需要分析方法的计算量有多大。

迭代法:

  • 构造一个收敛到解的迭代序列。产生的截断误差取决于迭代的步数,
  • 因而也没有截断误差的分析,但需要分析迭代格式收敛的条件。

先从直接法开始介绍。

在线性方程组$Ax=b$中,若$A$对角阵

\[\begin{pmatrix} a_{11} & 0 & \cdots & 0 \\ 0 & a_{12} & \cdots & 0 \\ \vdots & \vdots & & \vdots \\ 0 & 0 & \cdots & a_{nn} \\ \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{pmatrix} =\begin{pmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{pmatrix} \]

则可以直接得到解

\[x_i=\frac{b_i}{a_{ii}}, \quad i=1,2,\cdots,n \]

可以看到,计算得到解需要的计算量是: $n$个乘除法。

在线性方程组$Ax=b$中,若$A$为下三角阵

\[\begin{pmatrix} a_{11} & 0 & \cdots & 0 \\ a_{21} & a_{12} & \cdots & 0 \\ \vdots & \vdots & & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \\ \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{pmatrix} =\begin{pmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{pmatrix} \]

则可以顺序得到解

\[x_i=\frac{1}{a_{ii}}\left(b_i-\sum_{k=1}^{i-1}a_{ik}x_k\right), \quad i=1,2,\cdots,n \]

$x_i$需要运算量是: $i-1$个加减,$i$个乘除。因此,总的计算量是

\[\sum_{i=1}^n(2i-1)=n^2 \]

在线性方程组$Ax=b$中,若$A$上三角阵

\[\begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ 0 & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & & \vdots \\ 0 & 0 & \cdots & a_{nn} \\ \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{pmatrix} =\begin{pmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{pmatrix} \]

则可以倒序得到解

\[x_i=\frac{1}{a_{ii}}\left(b_i-\sum_{k={i+1}}^{n}a_{ik}x_k\right), \quad i=n,n-1,\cdots,1 \]

总计算量与下三角阵的情形一样,也是$n^2$次运算。

定义 1.
消元法就是把矩阵变为对角阵、上三角阵或下三角阵,然后进行求解的方法

Gauss 消元

Gauss消元法把矩阵

\[\begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1n} & | & b_1 \\ a_{21} & a_{22} & \cdots & a_{2n} & | & b_2 \\ \vdots & \vdots & \ddots & \vdots & | & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} & | & b_n \\ \end{pmatrix} \]

变为上三角阵

\[\begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1n} & | & b_1 \\ 0 & a_{22}^{(2)} & \cdots & a_{2n}^{(2)} & | & b_2^{(2)} \\ \vdots & \vdots & \ddots & \vdots & | & \vdots \\ 0 & 0 & \cdots & a_{nn}^{(n)} & | & b_n^{(n)} \\ \end{pmatrix} \]

第1步,第1行$\times\frac{-a_{i1}}{a_{11}}$加到第$i$行,$i=2,\cdots,n$,可以把第1列下三角部分变成$0$

\[\begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1n} & | & b_1 \\ a_{21} & a_{22} & \cdots & a_{2n} & | & b_2 \\ \vdots & \vdots & \ddots & \vdots & | & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} & | & b_n \\ \end{pmatrix} \rightarrow \begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1n} & | & b_1 \\ 0 & a_{22}^{(2)} & \cdots & a_{2n}^{(2)} & | & b_2^{(2)} \\ \vdots & \vdots & \ddots & \vdots & | & \vdots \\ 0 & a_{n2}^{(2)} & \cdots & a_{nn}^{(2)} & | & b_n^{(2)} \\ \end{pmatrix} \]

这一步的运算量$(n-1)\times(n-1+1+1)$次乘除,$(n-1)\times(n)$次加减, 共$(n-1)(2n+1)$次运算。

第2步,第2行$\times\frac{-a_{i2}^{(2)}}{a_{22}^{(2)}}$加到第$i$行,$i=3,\cdots,n$,可以把第2列下三角部分变成$0$

\[\begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1n} & | & b_1 \\ 0 & a_{22}^{(2)} & \cdots & a_{2n}^{(2)} & | & b_2^{(2)} \\ \vdots & \vdots & \ddots & \vdots & | & \vdots \\ 0 & a_{n2}^{(2)} & \cdots & a_{nn}^{(2)} & | & b_n^{(2)} \\ \end{pmatrix} \rightarrow \begin{pmatrix} a_{11} & a_{12} & \cdots & \cdots & a_{1n} & | & b_1 \\ 0 & a_{22}^{(2)} & \cdots & \cdots & a_{2n}^{(2)} & | & b_2^{(2)} \\ 0 & 0 & a_{33}^{(3)} & \cdots & a_{3n}^{(3)} & | & b_3^{(3)} \\ \vdots & \vdots & \vdots & \ddots & \vdots & | & \vdots \\ 0 & 0 & a_{n3}^{(3)} & \cdots & a_{nn}^{(3)} & | & b_n^{(3)} \\ \end{pmatrix} \]

这一步的运算量$(n-2)\times(n-2+1+1)$次乘除,$(n-2)\times(n-1)$次加减, 共$(n-2)(2n-1)$次运算。

在第$k-1$步后,

\[\begin{pmatrix} a_{11} & a_{12} &\cdots & a_{1k} & \cdots & a_{1n} & | & b_1 \\ & a_{22}^{(2)} & \cdots & a_{2k}^{(2)}& \cdots & a_{2n}^{(1)} & | & b_2^{(2)} \\ & & \ddots &\vdots & & \vdots & | & \vdots \\ & & & a_{kk}^{(k)} & \cdots & a_{kn}^{(k)} &|& b_k^{(k)} \\ & & & \vdots & \ddots & \vdots &|& \vdots \\ & & & a_{nk}^{(k)} & \cdots & a_{nn}^{(k)} &|& b_n^{(k)} \\ \end{pmatrix} \]

第k步,第k行$\times\frac{-a_{ik}^{(k)}}{a_{kk}^{(k)}}$加到第$i$行,$i=k+1,\cdots,n$,可以把第k列下三角部分变成$0$

这一步的运算量$(n-k)\times(2n-2k+3)$

经过$n-1$步后,可以得到上三角阵

\[\begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1n} & | & b_1 \\ & a_{22}^{(2)} & \cdots & a_{2n}^{(2)} & | & b_2^{(2)} \\ & & \ddots & \vdots & | & \vdots \\ & & & a_{nn}^{(n)} & | & b_n^{(n)} \\ \end{pmatrix} \]

加上一个解上三角矩阵的回代过程,即可得解$x_i$

整个消元过程的计算量为

\[\sum_{k=1}^{n-1}(n-k)(2n-2k+3)=\frac23n^3+\frac12n^2-\frac76n \]

加上反代过程的计算量为$n^2$,可以得到全部的运算量为

\[\frac23{n^3}+\frac32n^2-\frac76n=\frac23n^3+O(n^2) \]
  • 可以看到,整个Gauss消元过程都是对方程组作的初等变换,它不会改变方程组的根, 也不会改变$\det A$的值。因此,不存在截断误差。
  • Gauss消元法可以进行,需要保证
    \[a_{kk}^{(k)}\neq 0, \quad k=1,2,\cdots,n-1 \]

定理 1.
方程$Ax=b$可以使用Gauss消元法求解的充要条件是$A$的所有顺序主子式不为0,即

\[det\begin{pmatrix} a_{11} & \cdots & a_{1i} \\ \vdots & \ddots & \vdots \\ a_{i1} & \cdots & a_{ii} \end{pmatrix}\neq 0 , \forall i=1,2,\cdots,n \]

Gauss消元法的问题:

  • $\det A\neq 0$可以保证方程组存在唯一的解,但还不能保证Gauss消元法可以运行。
  • $|a_{kk}^{(k)}|$接近0时,舍入误差会被放大,从而有可能导致计算失败。

例 1. 解方程

\[\begin{cases} 10^{-8}x_1+x_2=1 \\ x_1+x_2=3 \end{cases} \]

. 方程有真解

\[\begin{cases} x_1=2.0000000200000002... \\ x_2=0.9999999799999997... \end{cases} \]

由Gauss消元

\[\begin{pmatrix} 10^{-8} & 1 & | & 1 \\ 1 & 1 & | & 3 \end{pmatrix} \rightarrow \begin{pmatrix} 10^{-8} & 1 & | & 1 \\ 0 & 1-10^8 & | & 3-10^8 \end{pmatrix} \]

单精度数值计算下,

\[\begin{aligned} 1-10^8=-9999999=-0.9999999\times 10^{8}\approx-1.000000\times 10^8 \\ 3-10^8=-9999997=-0.9999997\times 10^{8}\approx-1.000000\times 10^8 \end{aligned} \]

则有

\[\rightarrow \begin{pmatrix} 10^{-8} & 1 & | & 1 \\ 0 & -10^8 & | & -10^8 \end{pmatrix} \]

进一步得到

\[x_2=1, x_1=0 \]

$x_1$与真解的误差很大

列主元(Partial Pivot)

Gauss消元法的问题集中在$a_{kk}^{(k)}$上,

  • 若它为0,则算法不可用;
  • 若它太接近0,则可能导致不稳定: 舍入误差会导致结果误差太大。

在第$k-1$步后,

\[\begin{pmatrix} a_{11} & a_{12} &\cdots & a_{1k} & \cdots & a_{1n} & | & b_1 \\ 0 & a_{22}^{(2)} & \cdots & a_{2k}^{(2)}& \cdots & a_{2n}^{(1)} & | & b_2^{(2)} \\ \vdots & \vdots & \ddots &\vdots & & \vdots & | & \vdots \\ 0 & 0 & \cdots & a_{kk}^{(k)} & \cdots & a_{kn}^{(k)} &|& b_k^{(k)} \\ \vdots & \vdots & \cdots & \vdots & \ddots & \vdots &|& \vdots \\ 0 & 0 & \cdots & a_{nk}^{(k)} & \cdots & a_{nn}^{(k)} &|& b_n^{(k)} \\ \end{pmatrix} \]

$|a_{jk}^{(k)}|$$|a_{kk}^{(k)}|,|a_{k+1,k}^{(k)}|,\cdots,|a_{nk}^{(k)}|$中最大的,则交换$j$行与$k$行。 然后进行第k步消元。

例 2. 解方程

\[\begin{cases} 10^{-8}x_1+x_2=1 \\ x_1+x_2=3 \end{cases} \]

. 由列主元消元

\[\begin{aligned} \begin{pmatrix} 10^{-8} & 1 & | & 1 \\ 1 & 1 & | & 3 \end{pmatrix} & \rightarrow \begin{pmatrix} 1 & 1 & | & 3 \\ 10^{-8} & 1 & | & 1 \\ \end{pmatrix} \\ & \rightarrow \begin{pmatrix} 1 & 1 & | & 3 \\ 0 & 1-10^{-8} & | & 1-3\times 10^{-8} \\ \end{pmatrix} \end{aligned} \]

在单精度数值计算下,(同样存在大数吞小数)

\[1-10^{-8}=0.9999999\approx-1.000000 \]
\[1-3\times10^{-8}=0.9999997\approx-1.000000 \]

则有

\[\rightarrow \begin{pmatrix} 1 & 1 & | & 3 \\ 0 & 1 & | & 1 \end{pmatrix} \]

进一步得到

\[x_2=1, x_1=2 \]

与真解的误差在浮点数的精度范围内

例 3. 解方程

\[\begin{cases} x_1+10^8 x_2=10^8 \\ x_1+x_2=3 \end{cases} \]

. 由列主元消元

\[\begin{aligned} \begin{pmatrix} 1 & 10^8 & | & 10^8 \\ 1 & 1 & | & 3 \end{pmatrix} & \rightarrow \begin{pmatrix} 1 & 10^8 & | & 10^8 \\ 0 & 1-10^8 & | & 3-10^8 \\ \end{pmatrix} \\ \end{aligned} \]

在单精度数值计算下, 则有

\[\rightarrow \begin{pmatrix} 1 & 10^8 & | & 10^8 \\ 0 & -10^8 & | & -10^8 \end{pmatrix} \]

进一步得到$ x_2=1,\, x_1=0$

$x_1$的误差很大。

  • $\det A\neq 0$,则$|a_{kk}^{(k)}|,|a_{k+1,k}^{(k)}|,\cdots,|a_{nk}^{(k)}|$中 一定有一个非0元。(否则,$\det A=0$
  • 列主元消元的适用范围是$\det A\neq 0$,与线性方程组存在唯一解的条件是一致的。
  • 列主元在一定程序上克服了Gauss消元中的不稳定现象。
  • 相较于Gauss消元,列主元没有增加任何的运算量,但增加了$\displaystyle\sum_{i=1}^{n-1}(n-i)$次的比较操作。
  • 还存在一些情形,用列主元消元法仍然有较大的误差。

改进的方法可以有 Scaled Pivot全选主元消元

Scaled Pivot

在第$k-1$步后,

\[\begin{pmatrix} a_{11} & a_{12} &\cdots & a_{1k} & \cdots & a_{1n} & | & b_1 \\ 0 & a_{22}^{(2)} & \cdots & a_{2k}^{(2)}& \cdots & a_{2n}^{(1)} & | & b_2^{(2)} \\ \vdots & \vdots & \ddots &\vdots & & \vdots & | & \vdots \\ 0 & 0 & \cdots & a_{kk}^{(k)} & \cdots & a_{kn}^{(k)} &|& b_k^{(k)} \\ \vdots & \vdots & \cdots & \vdots & \ddots & \vdots &|& \vdots \\ 0 & 0 & \cdots & a_{nk}^{(k)} & \cdots & a_{nn}^{(k)} &|& b_n^{(k)} \\ \end{pmatrix} \]

$i=k,\cdots,n$中,比较

\[s_i=\frac{|a_{ik}^{(k)}|}{\displaystyle\max_{j=k,\cdots,n}|a_{ij}^{(k)}|} \]

$s_j$$s_k,s_{k+1},\cdots,s_n$中最大的,则交换$j$行与$k$行,然后进行消元。

整个算法,在第$k$步,增加了$n-k$个查找最大值的过程,以及$n-k$次除法。

因此,整个算法的运算量增加了$\frac{n(n-1)}2$次。

例 4. 解方程

\[\begin{cases} x_1+10^8 x_2=10^8 \\ x_1+x_2=3 \end{cases} \]

. 由Scaled Pivot

\[\begin{aligned} \begin{pmatrix} 1 & 10^8 & | & 10^8 \\ 1 & 1 & | & 3 \end{pmatrix} & \rightarrow \begin{pmatrix} 1 & 1 & | & 3 \\ 1 & 10^8 & | & 10^8 \\ \end{pmatrix} \\ & \rightarrow \begin{pmatrix} 1 & 1 & | & 3 \\ 1 & 10^8-1 & | & 10^8-3 \\ \end{pmatrix} \end{aligned} \]

在单精度数值计算下, 则有

\[\rightarrow \begin{pmatrix} 1 & 1 & | & 3 \\ 1 & 10^8 & | & 10^8 \\ \end{pmatrix} \rightarrow \begin{cases} x_1=2 \\ x_2=1 \end{cases} \]

全选主元 Full Pivot

在第$k-1$步后,

\[\begin{pmatrix} a_{11} & a_{12} &\cdots & a_{1k} & \cdots & a_{1n} & | & b_1 \\ 0 & a_{22}^{(2)} & \cdots & a_{2k}^{(2)}& \cdots & a_{2n}^{(2)} & | & b_2^{(2)} \\ \vdots & \vdots & \ddots &\vdots & & \vdots & | & \vdots \\ 0 & 0 & \cdots & a_{kk}^{(k)} & \cdots & a_{kn}^{(k)} &|& b_k^{(k)} \\ \vdots & \vdots & \cdots & \vdots & \ddots & \vdots &|& \vdots \\ 0 & 0 & \cdots & a_{nk}^{(k)} & \cdots & a_{nn}^{(k)} &|& b_n^{(k)} \\ \end{pmatrix} \]

$|a_{jl}^{(k)}|$$\{|a_{ij}^{(k)}|,i=k,\cdots,n$, $j=k,\cdots,n\}$中最大,则交换$k$行与$j$行,$k$列与$l$列。 然后进行第k步消元。

问题. 列交换并不是等价的变换,还有可能得到真正的解吗?

比较如下两个线性方程组

\[\begin{cases} a_{11}x_1+a_{12}x_2=b_1 \\ a_{21}x_1+a_{22}x_2=b_2 \end{cases} ,\quad \begin{cases} a_{12}y_1+a_{11}y_2=b_1 \\ a_{22}y_1+a_{21}y_2=b_2 \end{cases} \]

可以看出,

\[\begin{cases} x_1=y_2 \\ x_2=y_1 \end{cases} \]

这样,对增广矩阵的一次列交换,需要对结果做相应的一次行交换

全选主元消元法的运算量与Gauss消元列主元消元的运算量是一样的。

相对于Gauss消元,全选主元消元法增加了$\displaystyle\sum_{i=1}^{n-1}[(n-i+1)^2-1]$ 次的比较操作。

例 5. 解方程

\[\begin{cases} x_1+10^8 x_2=10^8 \\ x_1+x_2=3 \end{cases} \]

. 由全选主元消元

\[\begin{aligned} \begin{pmatrix} 1 & 10^8 & | & 10^8 \\ 1 & 1 & | & 3 \end{pmatrix} & \rightarrow \begin{pmatrix} 10^8 & 1 & | & 10^8 \\ 1 & 1 & | & 3 \end{pmatrix} \\ & \rightarrow \begin{pmatrix} 10^8 & 1 & | & 10^8 \\ 0 & 1-10^{-8} & | & 2 \end{pmatrix} \end{aligned} \]

在单精度数值计算下, 则有

\[\rightarrow \begin{pmatrix} 10^8 & 1 & | & 10^8 \\ 0 & 1 & | & 2 \end{pmatrix} \rightarrow \begin{cases} x_1=1 \\ x_2=2 \end{cases} \]

消元过程中,交换了1列与2列,因此,最后的解需要交换1行与2行。

得到解是

\[\begin{cases} x_1=2 \\ x_2=1 \end{cases} \]

Gauss-Jordan消元

消元过程中,把矩阵的上三角部分也变成0元,最后得到一个对角阵。

第k步: 第k行$\times\frac{-a_{ik}^{(k)}}{a_{kk}^{(k)}}$加到第$i$行,

$\quad\quad \quad \quad \quad i=1,\cdots,k-1,k+1,\cdots,n$

可以把第k列上三角部分与下三角部分变成$0$

  • 具体实现时,可以先做Gauss消元,变成上三角阵,然后从第$n$列开始,逐步将上三角部分也变成0。
  • $AX=B$形问题,在反代过程中,可以节省运算量。
  • $B=I$时,得到的$X=A^{-1}$

程序作业

编写列主元消元法,并且计算

\[\begin{pmatrix} 31 & -13 & 0 & 0 & 0 & -10 \\ -13 & 35 & -9 & 0 & -11 \\ & -9 & 31 & -10 \\ & & -10 & 79 & -30 & 0 & 0 & 0 & -9 \\ & & &-30 & 57 & -7 & 0 & -5 \\ & & & & -7 & 47 & -30 \\ & & & & & -30 & 41 \\ & & & & -5 & & & 27 & -2 \\ & & & -9 & & & &-2 & 29 \\ \end{pmatrix} x= \begin{pmatrix} -15 \\ 27 \\ -23 \\ 0 \\ -20 \\ 12 \\ -7 \\ 7 \\ 10 \end{pmatrix} \]

输出消元后的上三角阵,及解

输出:

消元后的矩阵为
31,-13,0,0,0,-10,0,0,0
0, *,
...

解为:
0.218304
1.29484
...

谢谢

vertical slide 2

目录

本节读完

例 6.

[#ex9-1-0].