线性方程组数值解法

§1 Gauss 消去法,直接三角分解解法 §2 Jacobi 迭代法,Gauss-Seidel 迭代法* §3 松弛迭代法 §4 极小化方法#

在自然科学和工程技术中很多问题的解决常常归结为求解线性方程组,例如求三次样条函数、最小二乘法拟合数据、有限差分或有限元求解偏微分方程边值问题等。由于克莱默法则求解线性方程组计算量大,效率低,因此,线性方程组的有效数值解法成为了非常重要且实用的工具。本章介绍求解各类线性方程组的常用数值方法,包括求解相容线性方程组的直接法和迭代法、求解矛盾方程组的广义逆矩阵法等。

Gauss 消去法

求解线性方程组的直接法指不计舍入误差的情况下,经过有限步算术运算,求得方程组精确解的方法。但实际计算中由于舍入误差的存在,这种方法也只能求得近似解。Gauss消去法就是求解线性方程组的一种直接法。它的主要思想是通过逐步消元(行初等变换),将原方程组化为等价的上三角方程组,然后回代求解。对于主元很小的情况,Gauss消去法可能产生很大的舍入误差,为避免这一缺点,可对其改进得到列主元法。

本节介绍Gauss消去法和列主元法。

设 \(\boldsymbol{A} = (a_{ij}) \in \mathbb{C}^{n \times n}\) , \(\boldsymbol{b} = [b_1, b_2, \cdots, b_n]^T \in \mathbb{C}^n\),求解线性方程组

\[\boldsymbol{Ax} = \boldsymbol{b} \tag{1}\]

如果 \(\boldsymbol{A}\) 是上三角矩阵或者下三角矩阵,我们都可以逐次回代求出所有的分量。

Gauss消去法的基本思想:将系数矩阵 \(\boldsymbol{A}\) 逐次消去成上三角矩阵,再逐次向后回代求出方程组所有的解。

高斯消去法中每步要用到主元素 \(a_{ii}\) 作为除数,当主元素绝对值很小时,计算机浮点数精度有限,就会有很大的误差。

为了解决这个问题,使用列主元素法

三角分解解法

三角分解解法是求解线性方程组的一种直接法,也是Gauss消去法的另一种表示形式。三角分解法的主要思想是通过矩阵的三角分解,将系数矩阵分解为一个下三角阵和一个上三角阵的乘积,然后分别用前代法和回代法求解相应的上三角方程组和下三角方程组。本节介绍基于不同三角分解的Doolittle法、适用于Hermite正定阵的Cholesky方法和改进的平方根法。

广义逆矩阵求解矛盾方程组

如果对于任意 \(\boldsymbol{x} \in \mathbb{C}^n\),都有 \(\boldsymbol{Ax} - \boldsymbol{b} \ne \boldsymbol{0}\) ,这时候 \(\boldsymbol{Ax} = \boldsymbol{b}\) 叫做矛盾方程组,对于这样的方程组,希望找到一个 \(\boldsymbol{x}_0\) 使得 \(\Vert \boldsymbol{Ax} - \boldsymbol{b} \Vert_2\) 最小,即

\[\Vert \boldsymbol{Ax}_0 - \boldsymbol{b} \Vert_2 = \min_{\boldsymbol{x} \in \mathbb{C}^n} \Vert \boldsymbol{Ax} - \boldsymbol{b} \Vert_2\]

这个问题就是线性方程组 \(\boldsymbol{Ax} = \boldsymbol{b}\) 的最小二乘问题,此解为矛盾方程组的最小二乘解。

矛盾方程组,无法采用 Gauss 消去法或三角分解解法求解,但是,可以根据广义逆矩阵的相关知识,表示出矛盾方程组的最小二乘解极小范数最小二乘解。本节主要介绍利用广义逆矩阵求解矛盾方程组的方法,以及列满秩最小二乘问题的求解。

迭代法的一般概念

与直接法不同,迭代法一般不能通过有限次的算术运算求得线性方程组的精确解,而是从初始向量出发,利用设计好的步骤逐步逼近线性方程组的精确解,直至达到精度要求。随着计算机的飞速发展和问题规模的扩大,迭代法已成为求解大型稀疏线性方程组最重要的一类方法。本节主要介绍迭代法的一般概念,包括迭代格式、收敛性及其判定、收敛速度、误差估计等。

一个线性方程组 \(\boldsymbol{Ax} = \boldsymbol{b}\) ,变成同解方程组 \(\boldsymbol{x} = \boldsymbol{Bx} + f\),这就是一个迭代格式

\[\boldsymbol{x}^{(k+1)} = \boldsymbol{Bx}^{(k)} + f\]

其中 \(\boldsymbol{B}\) 为迭代矩阵。如果迭代序列满足 \(\lim_{k \to \infty} \boldsymbol{x}^{(k)} = \boldsymbol{x}^*\),那么迭代法收敛。

迭代不一定总是收敛,为什么说一个方法是可以用的呢?算法背后的数学原理,迭代收敛的等价命题:

  • 迭代法 \(\boldsymbol{x} = \boldsymbol{Bx} + f\) 收敛
  • 迭代矩阵 \(\rho(\boldsymbol{B}) < 1\)
  • 至少存在一个矩阵范数 \(\Vert \centerdot \Vert\),使得 \(\Vert \boldsymbol{B} \Vert < 1\)

迭代法的构造思路

系数矩阵可以拆解为

\[\boldsymbol{A} = \boldsymbol{M} - \boldsymbol{N} \tag{1}\]

其中 \(\boldsymbol{M}\) 非奇异。则由 \(\boldsymbol{Ax} = \boldsymbol{b}\) 可得

\[\boldsymbol{x} = \boldsymbol{M}^{-1}\boldsymbol{Nx} + \boldsymbol{M}^{-1}\boldsymbol{b}\]

这就有了 \(\boldsymbol{x} = \boldsymbol{Bx} + f\) 的形式

Jacobi 迭代法和 Gauss-Seidel 迭代法

迭代法求解线性方程组的主要思想是,将原方程组转化为具有某种形式的同解方程组,然后据此构造迭代格式。本节介绍求解线性方程组的两种简单迭代法:J迭代法和G-S迭代法。包括迭代格式的构造、迭代法的分量形式和矩阵形式以及迭代法敛散性的判定。

非奇异线性方程组 \(\boldsymbol{Ax} = \boldsymbol{b} , a_{ii} \ne 0\) ,把 \(\boldsymbol{A}\) 拆成 \(\boldsymbol{A} = \boldsymbol{D} + \boldsymbol{L} + \boldsymbol{U}\) ,其中

\[\boldsymbol{D} = \text{diag}\{ a_{11}, a_{22}, \cdots ,a_{nn}\}\] \[\boldsymbol{L} = \left[ \begin{array}{} 0 \\ a_{21} & 0 \\ a_{31} & a_{32} & 0 \\ \vdots & \vdots & \ddots & \ddots\\ a_{n1} & a_{n2} & \cdots & a_{n,(n-1)} & 0 \end{array} \right]\] \[\boldsymbol{U} = \left[ \begin{array}{} 0 & a_{12} & a_{13} & \cdots & a_{1n}\\ & 0 & a_{23} & \cdots & a_{2n}\\ & & \ddots & \ddots & \vdots\\ & & & 0 & a_{(n-1),n}\\ & & & & 0 \end{array} \right]\]

令 \(\boldsymbol{M} = \boldsymbol{D}\), \(\boldsymbol{N} = -\boldsymbol{L} - \boldsymbol{U}\) 得到的迭代式是 Jacobi 迭代式,如果 \(\boldsymbol{M} = \boldsymbol{D} + \boldsymbol{L}\), \(\boldsymbol{N} = - \boldsymbol{U}\),得到的就是 Gauss-Seidel 迭代式

Jacobi 迭代法的迭代矩阵

\[\] \[\boldsymbol{B} = -\boldsymbol{(D+L)^{-1}U }\]

超松弛迭代法

Jacobi 迭代法和 Gauss-Seidel 迭代法的迭代格式简单,但收敛速度可能很慢。为了提高迭代法的收敛速度,对 Gauss-Seidel 迭代法加以改进,便得到超松弛迭代法。本节介绍求解线性方程组的松弛迭代法,包括迭代格式的构造、迭代法的分量形式和矩阵形式以及迭代法敛散性的判定。

极小化方法

求解线性方程组的极小化方法也属于迭代法。它首先借助变分法的思想,将线性方程组的求解转化为二次函数的极值问题,然后构造不同的"下山"算法逼近该极值问题的解。根据"下山"方向选取的不同,本节介绍最速下降法和共轭梯度法。其中,最速下降法是仅具有局部最优性的算法,可能收敛很慢,而共轭梯度法是具有某种全局最优性的算法,在不计舍入误差情况下,至多进行n步便得方程组的精确解。