第二章 线性方程组的直接方法(2)

上一讲我们介绍了Gauss消去法的局限性，为了避免其局限性，提高计算的数值稳定性，在消元过程中采用选择主元的方法。常采用的是列主元消去法和全主元消去法。

2.2 列主元Gauss消去法

给定线性方程组$Ax=b$ ，记$A^{(1)} = A,b^{(1)} = b$，列主元Gauss消去法的具体过程如下：

（1）首先在增广矩阵$B^{(1)} = (A^{(1)}, b^{(1)})$的第一列元素中，取：

在第一列中取绝对值最大的元素为主元素，然后将该行与第一行进行互换，即

$|a_{k1}^{(1)}| = max_{1\leq i\leq n}|a_{i1}^{(1)}| ， r_k \Leftrightarrow r_1$

然后进行第一步消元得增广矩阵$B^{(2)} = (A^{(2)}, b^{(2)})$。

（2）再在矩阵$B^{(2)} = (A^{(2)}, b^{(2)})$的第二列元素中，取：

$|a_{k2}^{(2)}| = max_{2\leq i\leq n}|a_{i2}^{(2)}| ， r_k \Leftrightarrow r_2$

然后进行第二步消元得增广矩阵$B^{(3)} = (A^{(3)}, b^{(3)})$。

（3）按此方法继续进行下去，经过$n-1$步选主元和消元运算，得到增广矩阵$B^{(n)} = (A^{(n)},b^{(n)})$。则方程组$A^{(n)}x = b^{(n)}$是与原方程等价的上三角形方程组，可进行回代求解。

列主元Gauss消去法可以很好地提高数值稳定性，并且每次选择一列中绝对值最大的数来作为主元，那么它可以执行完消元的充要条件是，

只要$|A| \not = 0$，列主元Gauss消去法就可以顺利进行。

例2：采用十进制四位浮点计算，分别用顺序Gauss消去法和列主元Gauss消去法求解线性方程组：

$\begin{cases} 0.012x_1+0.01x_2+0.167x_3 = 0.6781 \\ x_1 + 0.8334x_2 + 5.91x_3 = 12.1 \\ 3200x_1 + 1200x_2 + 4.2x_3 = 981 \end{cases}$

我们观察发现，这个线性方程组的数量级差距很大，因此在加减和乘除运算中就很容易出现问题，可能会产生“大数吃小数”的问题。

方程组具有四位有效数字的精确解为：$x_1^* = 17.46，x_2^* = -45.76, x_3^* = 5.546$

解：1. 用顺序Gauss消去法求解，消元过程为：

$\left[\begin{matrix} 0.0120 & 0.0100 & 0.167 & 0.6781 \\ 1.000 & 0.8334 & 5.910 & 12.10 \\ 3200 & 1200 & 4.200 & 981.0 \end{matrix}\right]$

$= \left[\begin{matrix} 0.0120 & 0.0100 & 0.167 & 0.6781 \\ 0 & 0.1000×10^{-3} & -8.010 & -44.41 \\ 0 & -1467 & -4453×10 & -1798×10^2 \end{matrix}\right]$

$= \left[\begin{matrix} 0.0120 & 0.0100 & 0.167 & 0.6781 \\ 0 & 0.1000×10^{-3} & -8.010 & -44.41 \\ 0 & 0 & -1175×10^5 & -6517×10^5 \end{matrix}\right]$

回代得$x_3 = 5.546, x_2 = 100.0, x_1 = -104.0$

我们前面说过，顺序Gauss消去法是直接方法，不存在截断误差，那么导致误差巨大的原因就来自于舍入误差，产生舍入误差大的原因在于主元非常小。

使用列主元Gauss消去法求解，消元过程为：

$\left[\begin{matrix} 0.0120 & 0.0100 & 0.167 & 0.6781 \\ 1.000 & 0.8334 & 5.910 & 12.10 \\ 3200 & 1200 & 4.200 & 981.0 \end{matrix}\right]$

选主元，将第一行和第三行互换，得：$\left[\begin{matrix} 3200 & 1200 & 4.200 & 981.0\\ 1.000 & 0.8334 & 5.910 & 12.10 \\ 0.0120 & 0.0100 & 0.167 & 0.6781 \end{matrix}\right]$

消元得：$\left[\begin{matrix} 3200 & 1200 & 4.200 & 981.0\\ 0 & 0.4584 & 5.909 & 11.79 \\ 0 & 0.55×10^{-2} & 0.167 & 0.6744 \end{matrix}\right]$

$\left[\begin{matrix} 3200 & 1200 & 4.200 & 981.0\\ 0 & 0.4584 & 5.909 & 11.79 \\ 0 & 0 & 0.0961 & 0.5329 \end{matrix}\right]$

回代得：$x_3 = 5.545, x_2 = -45.77, x_1 = 17.46$

因此我们看到，列主元Gauss消去法确实能够提高方法的数值稳定性。

小结：列主元Gauss消去法是在每一步消元前，在主元所在的一列选取绝对值最大的元素作为主元素。

而还有一种方法叫做全主元Gauss消去法，是在每一步消元前，在所有元素中选取绝对值最大的元素作为主元素。但由于运算量大增，实际应用中并不常用。

2.3 矩阵三角分解法

我们在消去时都是采用增广矩阵来进行消去，目的是为了保持方程的等价性。本节我们仅对系数矩阵来进行矩阵运算。每对矩阵进行一次消元，实际上就是实施了一次初等行变换。那么从左乘初等矩阵的角度来看是怎样的呢？

对矩阵$A^{(1)} = \left[\begin{matrix} a_{11}^{(1)} & a_{12}^{(1)} & a_{13}^{(1)}&... & a_{1n}^{(1)} \\ a_{21}^{(1)} & a_{22}^{(1)} & a_{23}^{(1)}& ... & a_{2n}^{(1)} \\a_{31}^{(1)} & a_{32}^{(1)} &a_{33}^{(1)} & ... & a_{3n}^{(1)}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ a_{n1}^{(1)} & a_{n2}^{(1)} & a_{n3}^{(1)} & ... & a_{nn}^{(1)}\end{matrix}\right]$

第一步：若$a_{11}^{(1)} \not = 0$，令$l_{i1} = \frac{a_{i1}^{(1)}}{a_{11}^{(1)}}, (i=2,3,...,n)$，记：

$L_1 = \left[\begin{matrix} 1 & & & & \\ -l_{21} & 1 & & & \\ -l_{31} & & 1 & & \\ \vdots & & & \ddots & \\ -l_{n1} & & & & 1\end{matrix}\right]$

像$L_1$这样的矩阵称为单位下三角矩阵，其主对角线上元素都是1，上三角部分都是0。接下来用$L_1$左乘$A^{(1)}$得：

$A^{(2)} = L_1A^{(1)} = \left[\begin{matrix} a_{11}^{(1)} & a_{12}^{(1)} & a_{13}^{(1)}&... & a_{1n}^{(1)} \\ 0 & a_{22}^{(2)} & a_{23}^{(2)}& ... & a_{2n}^{(2)} \\ 0 & a_{32}^{(2)} &a_{33}^{(2)} & ... & a_{3n}^{(2)}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\0 & a_{n2}^{(2)} & a_{n3}^{(2)} & ... & a_{nn}^{(2)}\end{matrix}\right]$

第二步：若$a_{22}^{(2)} \not = 0$，令$l_{i2} = \frac{a_{i2}^{(2)}}{a_{22}^{(2)}}, (i=3,4,...,n)$，记：

$L_2 = \left[\begin{matrix} 1 & & & & \\ & 1 & & & \\ & -l_{32} & 1 & & \\ & \vdots & & \ddots & \\ & -l_{n2} & & & 1\end{matrix}\right]$

矩阵$L_2$仍然是单位下三角矩阵，我们用$L_2$去左乘矩阵$A^{(2)}$，得：

$A^{(3)} = L_2A^{(2)} = \left[\begin{matrix} a_{11}^{(1)} & a_{12}^{(1)} & a_{13}^{(1)}&... & a_{1n}^{(1)} \\ 0 & a_{22}^{(2)} & a_{23}^{(2)}& ... & a_{2n}^{(2)} \\ 0 & 0 &a_{33}^{(2)} & ... & a_{3n}^{(2)}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\0 & 0 & a_{n3}^{(2)} & ... & a_{nn}^{(2)}\end{matrix}\right]$

如此进行下去，第$n-1$步得到：

$A^{(n)} = L_{n-1}A^{(n-1)} = \left[\begin{matrix} a_{11}^{(1)} & a_{12}^{(1)} & a_{13}^{(1)}&... & a_{1n}^{(1)} \\ & a_{22}^{(2)} & a_{23}^{(2)}& ... & a_{2n}^{(2)} \\ & &a_{33}^{(2)} & ... & a_{3n}^{(2)}\\ & & & \ddots & \vdots\\ & & & & a_{nn}^{(2)}\end{matrix}\right]$

整理整个过程得：

$A^{(n)} = L_{n-1}A^{(n-1)} = L_{n-1}L_{n-2}A^{(n-2)} = ... = L_{n-1}L_{n-2}...L_2L_1A^{(1)}$

其中每一个$L_k$都是一个单位下三角矩阵，而且都是可以求解出来的，说明它们都是可逆矩阵，所以有：

$A = A^{(1)} = L_{1}^{-1}L_{2}^{-1}...L_{n-1}^-1A^{(n)} = LU$

这就是$A$的LU三角分解。

$A = LU = \left[\begin{matrix} 1 & 0 & 0 & 0 \\ l_{21} & 1 & 0 & 0 \\ l_{31} & l_{32} & 1 & 0 \\ l_{41} & l_{42} & l_{43} & 1\end{matrix}\right]\left[\begin{matrix} u_{11} & u_{12} & u_{13} & u_{14} \\ 0 & u_{22} & u_{23} & u_{24} \\ 0 & 0 & u_{33} & u_{34} \\ 0 & 0 & 0 & u_{44}\end{matrix}\right]$

2.4 直接三角分解法

定理1：若$n$阶方阵$A$的各阶顺序主子式不为0，则存在唯一单位下三角矩阵$L$和上三角矩阵$U$使$A=LU$。

证明：只证唯一性，设有两种分解$A=LU=L'U'$

则有$L'^{-1}L = U'U^{-1} = I$（若下三角矩阵=上三角矩阵，则等号两边都为单位矩阵）

所以$L=L',U=U'$。

于是$Ax=b$可以写成$LUx=b$，令$Ux=y$，得：

$\begin{cases} Ly=b\\Ux=y \end{cases}$

大家知道，矩阵乘向量得到的是向量，那么$Ly=b$，$Ux=y$。那么原来一个线性方程组的求解，经过三角分解后变成两个三角形方程组的求解，而三角形线性方程组的求解，相对于普通线性方程组的求解要容易的很多。我们刚刚给出的三角分解和顺序Gauss消元法是一样的。

下面介绍三角分解的Doolittle分解方法，设：

$\left[\begin{matrix} a_{11} & a_{12} & ... & a_{1n} \\ a_{21} & a_{22} & ... & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & ... & a_{nn} \end{matrix}\right] = \left[\begin{matrix} 1 & & & &\\ l_{21} & 1 & & & \\ l_{31} & l_{32} & 1 & & \\ \vdots & \vdots & \vdots & \ddots & \\ l_{n1} & l_{n2} & l_{n3}& ... & 1\end{matrix}\right]\left[\begin{matrix} u_{11} & u_{12} & ...& u_{1n} \\ & u_{22} & ... & u_{2n} \\ & & \ddots & \vdots \\ & & & u_{nn}\end{matrix}\right]$

为了求得两个三角矩阵的分量元素，我们要利用两个三角矩阵的特殊形式和矩阵乘法的运算规则。

则$\begin{cases} u_{1j} = a_{1j}, j=1,2,...,n \\ l_{i1} = a_{i1}\div u_{11}, i=2,3,...,n \\ a_{kj} = l_{k1}u_{1j}+l_{k2}u_{2j}+...+l_{kk-1}u_{k-1j}+u_{kj} \\ u_{kj} = a_{kj}-\sum ^{k-1}_{m=1}l_{km}u_{mj}, j=k,k+1,...,n\\ l_{ik} = (a_{ik} -\sum ^{k-1}_{m=1}l_{im}u_{mk})\div u_{kk}, i=k+1,k+2,...,n \end{cases}$

求解顺序：先求$U$的第一行，再求$L$的第一列，再求$U$的第二行，再求$L$的第二列，以此类推。

由$\begin{cases} Ly=b\\Ux=y \end{cases}$，得

$\left[\begin{matrix} 1 & & & &\\ l_{21} & 1 & & & \\ l_{31} & l_{32} & 1 & & \\ \vdots & \vdots & \vdots & \ddots & \\ l_{n1} & l_{n2} & l_{n3}& ... & 1\end{matrix}\right] \left[\begin{matrix}y_1\\y_2\\\vdots\\y_n\end{matrix}\right] = \left[\begin{matrix}b_1\\b_2\\\vdots\\b_n\end{matrix}\right]$

$\left[\begin{matrix} u_{11} & u_{12} & ...& u_{1n} \\ & u_{22} & ... & u_{2n} \\ & & \ddots & \vdots \\ & & & u_{nn}\end{matrix}\right] \left[\begin{matrix}x_1\\x_2\\\vdots\\x_n\end{matrix}\right] = \left[\begin{matrix}y_1\\y_2\\\vdots\\y_n\end{matrix}\right]$

可得：

$\begin{cases} y_1=b_1 \\ y_{k} = b_{k} - \sum^{k-1}_{i=1} l_{ki}y_i, k=2,3,...,n \\ x_n = y_n \div u_{nn} \\ x_{i} = (y_i - \sum ^n_{j=i+1}u_{ij}x_j)\div u_{ii}, i=n-1,n-2,...,1 \end{cases}$

这就是求解方程组$Ax=b$的Doolittle三角分解法。

应用

例1：利用三角分解方法求解线性方程组

$\begin{cases} x_1+2x_2-3x_3 = 1 \\ 2x_1 - x_2 + 3x_3 = 5 \\ 3x_1 - 2x_2 + 2x_3 = 1 \end{cases}$

解：

(1) 将系数矩阵$A$书写出来为：

$A = \left[\begin{matrix} 1 & 2 & -3 \\ 2 & -1 & 3 \\ 3 & -2 & 2 \end{matrix}\right]$

(2) 将系数矩阵$A$分解为$LU$形式，分解后$U$的第一行和$A$的第一行一致，$L$的第一列和$A$的第一列一致，我们只需要计算其他位置的元素即可，有：

$A = \left[\begin{matrix} 1 & 2 & -3 \\ 2 & -1 & 3 \\ 3 & -2 & 2 \end{matrix}\right] \Rightarrow \left[\begin{matrix} 1 & 2 & -3 \\ 2 & -5 & 9 \\ 3 & \frac{8}{5} & -\frac{17}{5} \end{matrix}\right]$

所以，我们就可以写出$L$和$U$了

$A = \left[\begin{matrix} 1 & 2 & -3 \\ 2 & -1 & 3 \\ 3 & -2 & 2 \end{matrix}\right] = \left[\begin{matrix} 1 & & \\ 2 & 1 & \\ 3 & \frac{8}{5} & 1\end{matrix}\right]\left[\begin{matrix} 1 & 2 & -3 \\ & -5 & 9 \\ & & -\frac{17}{5} \end{matrix}\right] = LU$

(3) 接下来就是求解两个三角方程组了，先解$Ly=b$，有：

$\left[\begin{matrix} 1 & & \\ 2 & 1 & \\ 3 & \frac{8}{5} & 1\end{matrix}\right]\left[\begin{matrix} y_1 \\ y_2 \\ y_3 \end{matrix}\right] = \left[\begin{matrix} 1 \\ 5 \\ 1 \end{matrix}\right]$，得$\left[\begin{matrix} y_1 \\ y_2 \\ y_3 \end{matrix}\right] = \left[\begin{matrix} 1 \\ 3 \\ -\frac{34}{5} \end{matrix}\right]$

(4) 再解$Ux=y$，有：

$\left[\begin{matrix} 1 & 2 & -3 \\ & -5 & 9 \\ & & -\frac{17}{5} \end{matrix}\right]\left[\begin{matrix} x_1 \\ x_2 \\ x_3 \end{matrix}\right] = \left[\begin{matrix} 1 \\ 3 \\ -\frac{34}{5} \end{matrix}\right]$，得$\left[\begin{matrix} x_1 \\ x_2 \\ x_3 \end{matrix}\right] = \left[\begin{matrix} 1 \\ 3 \\ 2 \end{matrix}\right]$

解线性方程组$Ax=b$的Doolittle三角分解法的计算量约为$\frac{1}{3}n^3$，与Gauss消去法基本相同。其优点在于求一系列同系数的线性方程组$Ax=b_k(k=1,2,..,m)$时，可大大节省运算量。我们做三角分解仅仅是对系数矩阵$A$进行分解，而$b_k$的变化体现在后面求解两个三角方程组，所以我们说，如果用Doolittle三角分解法在求同系数三角分解时，只需要做一次分解。