概论 误差

误差的来源

模型误差

• 用计算机解决科学计算问题首先要建立数学模型，它是对被描述的实际问题进行抽象、简化而得到的，因而是近似的。我们把数学模型与实际问题之间出现的这种误差称为模型误差。
• 由于这种误差难以用数量表示，通常都假定数学模型是合理的，这种误差可忽略不计

观测误差

在数学模型中往往还有一些根据观测得到的物理量，如温度、长度、电压等，这些参量显然也包含误差。这种由观测产生的误差称为观测误差。

截断误差

当数学模型不能得到精确解时，通常要用数值方法求它的近似解，其近似解与精确解之间的误差称为截断误差或方法误差

• 例如，当函数 $f(x)$ 用 Taylor 多项式 $$P_{n}(x)=f(0)+\frac{f^{\prime}(0)}{1!} x+\frac{f^{\prime \prime}(0)}{2!} x^{2}+\cdots+\frac{f^{(n)}(0)}{n!} x^{n}$$ 近似代替时，数值方法的截断误差是 $$R_{n}(x)=f(x)-P_{n}(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!} x^{n+1},\quad \xi \text{在} x \text{与} 0 \text{之间}$$

舍入误差

有了求解数学问题的计算公式以后，用计算机进行数值计算时，由于计算机的字长有限，原始数据在计算机上表示会产生误差，计算过程又可能产生新的误差，这种误差称为舍入误差，例如，用 $3.14159$ 近似代替 $\pi$，产生的误差就是舍入误差

$$R=\pi-3.14159=0.0000026 \cdots$$

误差的基本概念

绝对误差

设 $x$ 为准确值，$x^{*}$ 为 $x$ 的一个近似值，称 $e^{*}=x^{*}-x$ 为近似值的绝对误差，简称误差

• 误差 $e^{*}$ 可正可负，当绝对误差为正时近似值偏大，叫做强近似值；当绝对误差为负时近似值偏小，叫做弱近似值

通常，我们不能算出准确值 $x$，也不能算出误差 $e^{*}$ 的准确值，只能根据测量工具或计算情况估计出误差的绝对值不超过某正数 $\varepsilon^{*}$，也就是误差绝对值的一个上界 $\varepsilon^{*}$，叫做近似值的误差限，它总是正数

对于一般情形，$\left|x^{*}-x\right| \leqslant \varepsilon^{*}$，即 $x^{*}-\varepsilon^{*} \leqslant x \leqslant x^{*}+\varepsilon^{*}$，这个不等式有时也表示为 $x=x^{*} \pm \varepsilon^{*}$

相对误差

除考虑误差的大小外，还应考虑准确值 $x$ 本身的大小，近似值的误差 $e^{*}$ 与准确值 $x$ 的比值

$$\frac{e^{*}}{x}=\frac{x^{*}-x}{x}$$

称为近似值 $x^{*}$ 的相对误差，记作 $e_{r}^{*}$

在实际计算中，由于真值 $x$ 总是不知道的，通常取

$$e_{\mathrm{r}}^{*}=\frac{e^{*}}{x^{*}}=\frac{x^{*}-x}{x^{*}}$$

作为 $x^{*}$ 的相对误差，条件是 $e_{r}^{*}=\frac{e^{*}}{x^{*}}$ 较小，此时

$$\frac{e^{*}}{x}-\frac{e^{*}}{x^{*}}=\frac{e^{*}\left(x^{*}-x\right)}{x^{*} x}=\frac{\left(e^{*}\right)^{2}}{x^{*}\left(x^{*}-e^{*}\right)}=\frac{\left(e^{*} / x^{*}\right)^{2}}{1-\left(e^{*} / x^{*}\right)}$$

是 $e_{r}^{*}$ 的二次方项级，故可忽略不计

相对误差也可正可负，它的绝对值上界叫做相对误差限，记作 $\varepsilon_{r}^{*},\varepsilon_{r}^{*}=\frac{\varepsilon^{*}}{\left|x^{*}\right|}$

有效数字

定义

若近似值 $x^{*}$ 的误差限是某一位的半个单位，该位到 $x^{*}$ 的第一位非零数字共有 $n$ 位，就说 $x^{*}$ 有 $n$ 位有效数字

• 如取 $x^{*}=3.14$ 作 $\pi$ 的近似值，$x^{*}$ 就有三位有效数字；取 $x^{*}=3.1416$ 作 $\pi$ 的近似值，$x^{*}$ 就有五位有效数字

标准形式

$x^{*}$ 有 $n$ 位有效数字可写成标准形式

$$x^{*}= \pm 10^{m} \times\left(a_{1}+a_{2} \times 10^{-1}+\cdots+a_{n} \times 10^{-(n-1)}\right)$$

其中，$a_{1}$ 是 $1$ 到 $9$ 中的一个数字；$a_{2},\cdots,a_{n}$ 是 $0$ 到 $9$ 中的一个数字；$m$ 为整数，且

$$\left|x-x^{*}\right| \leqslant \frac{1}{2} \times 10^{m-n+1}$$

和误差限的关系

以上述形式表示的有 $n$ 位有效数字的近似数 $x^{*}$，其绝对误差限为

$$\varepsilon^{*}=\frac{1}{2} \times 10^{m-n+1}$$

相对误差限为

$$\varepsilon_{r}^{*} \leqslant \frac{1}{2 a_{1}} \times 10^{-(n-1)}$$

反之，若 $x^{*}$ 的相对误差限

$$\varepsilon_{\mathrm{r}}^{*} \leqslant \frac{1}{2\left(a_{1}+1\right)} \times 10^{-(n-1)}$$

则 $x^{*}$ 至少具有 $n$ 位有效数字

数值运算的误差估计

两个近似数 $x_{1}^{*}$ 与 $x_{2}^{*}$，其误差限分别为 $\varepsilon\left(x_{1}^{*}\right)$ 及 $\varepsilon\left(x_{2}^{*}\right)$，它们进行加、减、乘、除运算得到的误差限分别为

$$\varepsilon\left(x_{1}^{*} \pm x_{2}^{*}\right)=\varepsilon\left(x_{1}^{*}\right)+\varepsilon\left(x_{2}^{*}\right)$$ $$\varepsilon\left(x_{1}^{*} x_{2}^{*}\right) \approx\left|x_{1}^{*}\right| \varepsilon\left(x_{2}^{*}\right)+\left|x_{2}^{*}\right| \varepsilon\left(x_{1}^{*}\right)$$ $$\varepsilon\left(\frac{x_{1}^{*}}{x_{2}^{*}}\right) \approx \frac{\left|x_{1}^{*}\right| \varepsilon\left(x_{2}^{*}\right)+\left|x_{2}^{*}\right| \varepsilon\left(x_{1}^{*}\right)}{\left|x_{2}^{*}\right|^{2}}\left(x_{2}^{*} \neq 0\right)$$

更一般的情况是，当自变量有误差时，计算函数值也会产生误差，其误差限可利用函数的 Taylor 展开式进行估计。
当 $f$ 为多元函数时计算 $A=f\left(x_{1},x_{2},\cdots,x_{n}\right)$，如果 $x_{1},x_{2},\cdots,x_{n}$ 的近似值为 $x_{1}^{*},x_{2}^{*} \cdots,x_{n}^{*}$，则 $A$ 的近似值为 $A^{*}=f\left(x_{1}^{*},x_{2}^{*},\cdots,x_{n}^{*}\right)$，由 Taylor 展开，得 $A^{*}$ 的相对误差限为

$$\varepsilon_{\mathrm{r}}^{*}=\varepsilon_{\mathrm{r}}\left(A^{*}\right)=\frac{\varepsilon\left(A^{*}\right)}{\left|A^{*}\right|} \approx \sum_{k=1}^{k}\left|\left(\frac{\partial f}{\partial x_{k}}\right)^{*}\right| \frac{\varepsilon\left(x_{k}^{*}\right)}{\left|A^{*}\right|}$$

数值运算中误差分析的方法与原则

运算过程舍入误差不增长的计算公式是数值稳定的，否则是不稳定的。

• 要避免除数绝对值远远小于被除数绝对值的除法
• 要避免两相近数相减
• 要防止大数“吃掉”小数
• 注意简化计算步骤，减少运算次数

第六章 方程求根

基本概念

数学、物理中的许多问题常常归结为求解函数方程 $f(x)=0$，这里，$f(x)$ 可以是代数多项式，也可以是超越函数。方程 $f(x)=0$ 的解 $x^{}$ 称为它的根，或称为 $f(x)$ 的*零点。

设函数 $f(x)$ 在 $[a,b]$ 上连续，且 $f(a) f(b)<0$，根据连续函数的性质可知方程 $f(x)=0$ 在区间 $(a,b)$ 内一定有实根，这时称 $[a,b]$ 为方程 $f(x)=0$ 的有根区间。

逐步搜索法

为明确起见，不妨假定 $f(a)<0,f(b)>0$。从有根区间 $[a,b]$ 的左端点 $x_{0}=a$ 出发，按某个预定的步长 $h$ （譬如取 $h=\frac{b-a}{N},N$ 为正整数）一步一步地向右跨，每跨一步进行一次根的“搜索”，即检查节点 $x_{k}=a+k h$ 上的函数值 $f\left(x_{k}\right)$ 的符号，一旦发现节点 $x_{k}$ 与端点 $a$ 的函数值异号，即 $f\left(x_{k}\right)>0^{\mathbb{1}}$，则可以确定一个缩小了的有根区间 $\left[x_{k-1},x_{k}\right]$，其宽度等于预定的步长 $h$

二分搜索法

再考察有根区间 $[a,b]$，取中点 $x_{0}=(a+b) / 2$ 将它分为两半，然后进行根的搜索，即检查 $f\left(x_{0}\right)$ 与 $f(a)$ 是否同号，如果确系同号，说明所求的根 $x^{*}$ 在 $x_{0}$ 的右侧，这时令 $a_{1}=x_{0},b_{1}=b$；否则 $x^{*}$ 必在 $x_{0}$ 的左侧，这时令 $a_{1}=a,b_{1}=x_{0}$。不管出现哪一种情况，新的有根区间 $\left[a_{1},b_{1}\right]$ 的长度仅为 $[a,b]$ 的一半。

对压缩了的有根区间 $\left[a_{1},b_{1}\right]$ 又可施行同样的手续，即用中点 $x_{1}=\left(a_{1}+b_{1}\right) / 2$ 将区间 $\left[a_{1},b_{1}\right]$ 再分为两半，然后通过根的搜索判定所求的根在 $x_{1}$ 的哪一侧，从而又确定一个新的有根区间 $\left[a_{2},b_{2}\right]$，其长度是 $\left[a_{1},b_{1}\right]$ 的一半。

如此反复二分下去，即可得出一系列有根区间

$$[a,b] \supseteq\left[a_{1},b_{1}\right] \supseteq\left[a_{2},b_{2}\right] \supseteq \cdots \supseteq\left[a_{k},b_{k}\right] \supseteq \cdots,$$

其中每个区间都是前一个区间的一半，因此 $\left[a_{k},b_{k}\right]$ 的长度 $b_{k}-a_{k}=(b-a) / 2^{k}$ 当 $k \rightarrow \infty$ 时趋于零。就是说，如果二分过程无限地继续下去，这些区间最终必将收缩于一点 $x^{*}$，该点显然就是所求的根。

迭代法

考察下列形式的方程：

$$x=\varphi(x)$$

这种方程是隐式的，因而不能直接得出它的根。但如果给出根的某个猜测值 $x_{0}$，将它代入上式的右端，即可求得 $x_{1}=\varphi\left(x_{0}\right)$。然后，又可取 $x_{1}$ 作为猜测值，进一步得到 $x_{2}=\varphi\left(x_{1}\right)$。如此反复迭代。如果按公式

$$x_{k+1}=\varphi\left(x_{k}\right),\quad k=0,1,2,\cdots$$

确定的数列 $\left\{x_{k}\right\}$ 有极限 $x^{*}=\lim _{k \rightarrow \infty} x_{k}$，则称迭代过程式收敛，这时极限值 $x^{*}$ 显然就是方程的根。

一般情况的收敛性

考察一般情形：设方程 $x=\varphi(x)$ 在区间 $(a,b)$ 内有根 $x^{*}$，则保证迭代过程 $x_{k+1}=\varphi\left(x_{k}\right)$ 收敛的条件可表述为：存在定数 $0<L<1$，使对于任意 $x \in[a,b]$ 成立

$$\left|\varphi^{\prime}(x)\right| \leqslant L$$

事实上，由微分中值定理

$$x_{k+1}-x^{*}=\varphi\left(x_{k}\right)-\varphi\left(x^{*}\right)=\varphi^{\prime}(\xi)\left(x_{k}-x^{*}\right),$$

其中 $\xi$ 是 $x^{*}$ 与 $x_{k}$ 之间某一点，当 $x_{k} \in[a,b]$ 时 $\xi \in[a,b]$。因此利用条件可以断定 $\left|x_{k+1}-x^{*}\right| \leqslant L\left|x_{k}-x^{*}\right|$。据此反复递推，有

$$\left|x_{k}-x^{*}\right| \leqslant L^{k}\left|x_{0}-x^{*}\right|$$

故当 $k \rightarrow \infty$ 时，迭代值 $x_{k}$ 将收敛到所求的根 $x^{*}$。

在上述论证过程中，应当保证一切迭代值 $x_{k}$ 全落在区间 $(a,b)$ 内，为此要求，对于任意 $x \in[a,b]$，总有 $\varphi(x) \in[a,b]$。

一般情况的误差估计

假定函数 $\varphi(x)$ 满足下列两项条件：

• 对于任意 $x \in[a,b]$，有 $$a \leqslant \varphi(x) \leqslant b$$
• 存在正数 $L<1$，使对于任意 $x \in[a,b]$，有 $$\left|\varphi^{\prime}(x)\right| \leqslant L<1$$ 则迭代过程 $x_{k+1}=\varphi\left(x_{k}\right)$ 对于任意初值 $x_{0} \in[a,b]$ 均收敛于方程 $x=\varphi(x)$ 的根 $x^{*}$，且有如下的误差估计式： $$\left|x_{k}-x^{*}\right| \leqslant \frac{L^{k}}{1-L}\left|x_{1}-x_{0}\right|$$

迭代过程是个极限过程，在用迭代法进行实际计算时，必须按精度要求控制送代次数。误差估计式原则上可用来确定迭代次数于实际应用。

$$\left|x^{*}-x_{k}\right| \leqslant \frac{1}{1-L}\left|x_{k+1}-x_{k}\right|$$

由此可见，只要相邻两次计算结果的偏差 $\left|x_{k+1}-x_{k}\right|$ 足够小，即可保证近似值 $x_{k}$ 具有足够的精度，因此可以通过检查 $\left|x_{k+1}-x_{k}\right|$ 来判断迭代过程应否终止。

误差估计（更强的条件）

假定函数 $\varphi(x)$ 满足下列两项条件：

• 对于任意 $x \in[a,b]$，有 $$a \leqslant \varphi(x) \leqslant b$$
• 存在正数 $L<1$，对任意 $x,y \in[a,b]$，有 $$\left|\varphi(x)-\varphi(y)\right| \leqslant L |x-y|$$ 则 $\varphi(x)$ 在 $[a,b]$ 存在唯一的不动点 $x^{*}$，这时称 $L$ 为收缩率。

计算步骤

下面列出迭代法的计算步骤：

1. 准备：提供迭代初值 $x_{0}$。
2. 迭代：计算迭代值 $x_{1}=\varphi\left(x_{0}\right)$。
3. 控制：检查 $\left|x_{1}-x_{0}\right|$：
   • 若 $\left|x_{1}-x_{0}\right|>\varepsilon$（$\varepsilon$ 为预先指定的精度），则以 $x_{1}$ 替换 $x_{0}$ 转步 $2$ 继续迭代
   • 当 $\left|x_{1}-x_{0}\right| \leqslant \varepsilon$ 时终止计算，取 $x_{1}$ 作为所求的结果。

在实际应用迭代法时，通常在所求的根 $x^{*}$ 的邻近进行考察，而研究所谓局部收敛性。

局部收敛性

• 若存在 $x^{*}$ 的某个邻域 $R:\left|x-x^{*}\right| \leqslant \delta$，使迭代过程 $x_{k+1}=\varphi\left(x_{k}\right)$ 对于任意初值 $x_{0} \in R$ 均收敛，则称迭代过程 $x_{k+1}=\varphi\left(x_{k}\right)$ 在根 $x^{*}$ 邻近具有局部收敛性。否则，若收敛范围为函数的定义域，称该迭代法具有全局收敛性。
• 设 $x^{*}$ 为方程 $x=\varphi(x)$ 的根，$\varphi^{\prime}(x)$ 在 $x^{*}$ 的邻近连续且 $\left|\varphi^{\prime}\left(x^{*}\right)\right|<1$，则迭代过程 $x_{k+1}=\varphi\left(x_{k}\right)$ 在 $x^{*}$ 邻近具有局部收敛性
• 若 $\varphi(x)$ 是 $[a,b]$ 上的连续函数，且 $\left|\varphi^{\prime}(x)\right|>1 \quad \forall x \in[a,b]$，则迭代 $x_{k+1}=\varphi\left(x_{k}\right)$ 在区间 $[a,b]$ 上不收敛。
• 设 $\varphi(x)$ 是 $[a,b]$ 上的连续函数，若 $x^{*}$ 是不动点方程 $x=\varphi(x)$ 的不动点，且 $\left|\varphi^{\prime}\left(x^{*}\right)\right|<1$，则 $\varphi(x)$ 在 $x^{*}$ 的邻域 $U_{\delta}\left(x^{*}\right)$ 内存在唯一的不动点 $x^{*},\forall x_{0} \in U_{\delta}\left(x^{*}\right) \quad x_{k+1}=\varphi\left(x_{k}\right) \rightarrow x^{*}$

收敛速度

对于一种迭代过程，为了保证它是有效的，需要肯定它的收敛性，同时考察它的收敛速度。所谓迭代过程的收敛速度，是指在接近收敛的过程中迭代误差的下降速度。

设迭代过程 $x_{k+1}=\varphi\left(x_{k}\right)$ 收敛于方程 $x=\varphi(x)$ 的根 $x^{*}$，如果迭代误差 $e_{k}=x_{k}-x^{*}$ 当 $k \rightarrow \infty$ 时成立下列渐近关系式：

$$\frac{e_{k+1}}{e_{k}^{p}} \rightarrow C \quad(C \neq 0 \mathrm{为常数})$$

则称该迭代过程是 p 阶收敛的

• 特别地，$p=1$ 时称为线性收敛，$p>1$ 时称为超线性收敛，$p=2$ 时称为平方收敛

对于迭代过程 $x_{k+1}=\varphi\left(x_{k}\right)$，如果 $\varphi^{(p)}(x)$ 在所求根 $x^{*}$ 的邻近连续，并且

$$\left\{\begin{array}{l} \varphi^{\prime}\left(x^{*}\right)=\varphi^{\prime \prime}\left(x^{*}\right)=\cdots=\varphi^{(p-1)}\left(x^{*}\right)=0,\\ \varphi^{(p)}\left(x^{*}\right) \neq 0, \end{array}\right.$$

则该迭代过程在点 $x^{*}$ 邻近是 $p$ 阶收敛的，且满足

$$\lim _{k \rightarrow \infty} \frac{\left(x_{k+1}-x^{*}\right)}{\left(x_{k}-x^{*}\right)^{p}}=\frac{\varphi^{(p)}\left(x^{*}\right)}{p!}$$

• 由上述定理可知，迭代过程的收敛速度依赖于迭代函数 $\varphi(x)$ 的选取.如果当 $x \in[a,b]$ 时 $\varphi^{\prime}(x) \neq 0$，则该迭代过程只可能是线性收敛的。
• Newton 法在根 $x_{*}$ 的临近是平方收敛的

牛顿迭代法

迭代函数：

$$\varphi(x)=x-\frac{f(x)}{f^{\prime}(x)},$$

Newton 公式

$$x\_{k+1}=x\_{k}-\frac{f\left(x\_{k}\right)}{f^{\prime}\left(x\_{k}\right)}$$

Newton 法有明显的几何解释。方程 $f(x)=0$ 的根 $x^{*}$ 可解释为曲线 $y=f(x)$ 与 $x$ 轴的交点的横坐标。设 $x_{k}$ 是根 $x^{*}$ 的某个近似值，过曲线 $y=f(x)$ 上横坐标为 $x_{k}$ 的点 $P_{k}$ 引切线，并将该切线与 $x$ 轴的交点的横坐标 $x_{k+1}$ 作为 $x^{*}$ 的新的近似值。注意到切线方程为

$$y=f\left(x_{k}\right)+f^{\prime}\left(x_{k}\right)\left(x-x_{k}\right)$$

从而就是 $Newton$ 公式的计算结果。由于这种几何背景，Newton 法亦称切线法

弦截法

设 $x_{k},x_{k-1}$ 是 $f(x)=0$ 的近似根，利用 $f\left(x_{k}\right),f\left(x_{k-1}\right)$ 构造一次插值多项式 $P_{1}(x)$，并用 $P_{1}(x)=0$ 的根作为 $f(x)=0$ 的新的近似根 $x_{k+1}$。由于

$$P_{1}(x) =f\left(x_{k}\right)+\frac{f\left(x_{k}\right)-f\left(x_{k-1}\right)}{x_{k}-x_{k-1}}\left(x-x_{k}\right)$$

因此有

$$x_{k+1} =x_{k}-\frac{f\left(x_{k}\right)}{f\left(x_{k}\right)-f\left(x_{k-1}\right)}\left(x_{k}-x_{k-1}\right)$$

这样导出的迭代公式可以看作 Newton 公式 $x_{k+1}=x_{k}-\frac{f\left(x_{k}\right)}{f^{\prime}\left(x_{k}\right)}$ 中的导数 $f^{\prime}(x)$ 用差商 $\frac{f\left(x_{k}\right)-f\left(x_{k-1}\right)}{x_{k}-x_{k-1}}$ 取代的结果。

• 切线法在计算 $x_{k+1}$ 时只用到前一步的值 $x_{k}$，而弦截法在求 $x_{k+1}$ 时要用到前面两步的结果 $x_{k},x_{k-1}$，因此使用这种方法必须先给出两个开始值 $x_{0},x_{1}$
• 收敛速度：$p=\frac{1+\sqrt{5}}{2} \approx 1.618$，超线性收敛

第七章 解线性方程组的直接方法

高斯（Gauss）消去法

高斯消去过程就是用行的初等变换将原方程组系数矩阵化为简单形式，从而将求解原方程组转化为求解简单方程组的问题。

基本过程

设有线性方程组

$$\left\{\begin{array}{c} a_{11} x_{1}+a_{12} x_{2}+\cdots+a_{1 n} x_{n}=b_{1},\\ a_{21} x_{1}+a_{22} x_{2}+\cdots+a_{2 n} x_{n}=b_{2},\\ \vdots \\ a_{n 1} x_{1}+a_{n 2} x_{2}+\cdots+a_{n n} x_{n}=b_{n}, \end{array}\right.$$

或写成矩阵形式 $\boldsymbol{A x}=\boldsymbol{b}$，其中

$$\boldsymbol{A}=\left(\begin{array}{cccc} a_{11} & a_{12} & \cdots & a_{1 n} \\ a_{21} & a_{22} & \cdots & a_{2 n} \\ \vdots & \vdots & & \vdots \\ a_{n 1} & a_{n 2} & \cdots & a_{n n} \end{array}\right),\quad \boldsymbol{x}=\left(\begin{array}{c} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \end{array}\right),\quad \boldsymbol{b}=\left(\begin{array}{c} b_{1} \\ b_{2} \\ \vdots \\ b_{n} \end{array}\right)$$

其中 $\boldsymbol{A}$ 为非奇异矩阵。

将方程记作 $\boldsymbol{A}^{(1)} \boldsymbol{x}=\boldsymbol{b}^{(1)}$，其中 $\boldsymbol{A}^{(1)}=\left(a_{i j}^{(1)}\right)=\left(a_{i j}\right),\quad \boldsymbol{b}^{(1)}=\boldsymbol{b}$

第一次消元：
设 $a_{11}^{(1)} \neq 0$，首先对行计算乘数 $m_{i 1}=a_{i 1}^{(1)} / a_{11}^{(1)}(i=2,3,\cdots,n)$，用 $-m_{i 1}$ 乘第 $1$ 个方程，加到第 $i(i=2,3,\cdots,n)$ 个方程上，消去第 $2$ 个方程直到第 $n$ 个方程中的未知数 $x_{1}$，得与原方程组等价的方程组

$$\left(\begin{array}{ccccc} a_{11}^{(1)} & a_{12}^{(1)} & \cdots & a_{1 n}^{(1)} \\ 0 & a_{22}^{(2)} & \cdots & a_{2 n}^{(2)} \\ \vdots & \vdots & & \vdots \\ 0 & a_{n 2}^{(2)} & \cdots & a_{n n}^{(2)} \end{array}\right)\left(\begin{array}{c} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \end{array}\right)=\left(\begin{array}{c} b_{1}^{(1)} \\ b_{2}^{(2)} \\ \vdots \\ b_{n}^{(2)} \end{array}\right)$$

简记作 $\boldsymbol{A}^{(2)} \boldsymbol{x}=\boldsymbol{b}^{(2)}$，其中

$$a_{i j}^{(2)}=a_{i j}^{(1)}-m_{i 1} a_{1 j}^{(1)},\quad b_{i}^{(2)}=b_{i}^{(1)}-m_{i 1} b_{1}^{(1)} \quad(i,j=2,3,\cdots,n)$$

一般第 $k(1 \leqslant k \leqslant n-1)$ 次消元：设第 $k-1$ 步计算已经完成，即已计算好与原方程组等价的方程组

$$\boldsymbol{A}^{(k)} \boldsymbol{x}=\boldsymbol{b}^{(k)}$$

且已消去未知数 $x_{1},x_{2},\cdots,x_{k-1}$，其中 $\boldsymbol{A}^{(k)}$ 具有如下形式：

$$\boldsymbol{A}^{(k)}=\left(\begin{array}{cccccc} a_{11}^{(1)} & a_{12}^{(1)} & \cdots & \cdots & \cdots & a_{1 n}^{(1)} \\ & a_{22}^{(2)} & \cdots & \cdots & \cdots & a_{2 n}^{(2)} \\ & & \ddots & & & \vdots \\ & & & a_{k k}^{(k)} & \cdots & a_{k n}^{(k)} \\ & & & \vdots & & \vdots \\ & & & a_{n k}^{(k)} & \cdots & a_{n n}^{(k)} \end{array}\right)$$

设 $a_{k k}^{(k)} \neq 0$，计算乘数 $m_{i k}=a_{i k}^{(k)} / a_{k k}^{(k)} \quad(i=k+1,\cdots,n)$，用 $-m_{i k}$ 乘上一次消去完的第 $k$ 个方程加上第 $i(i=k+1,\cdots,n)$ 个方程，消去第 $k+1$ 个方程直到第 $n$ 个方程的未知数 $x_{k}$，得到与原方程组等价的方程组 $\boldsymbol{A}^{(k+1)} \boldsymbol{x}=\boldsymbol{b}^{(k+1)}$
$\boldsymbol{A}^{(k+1)}$ 元素的计算公式为

$$\left\{\begin{array}{ll} a_{i j}^{(k+1)}=a_{i j}^{(k)}-m_{i k} a_{k i}^{(k)} & (i,j=k+1,\cdots,n),\\ b_{i}^{(k+1)}=b_{i}^{(k)}-m_{i k} b_{k}^{(k)} & (i=k+1,\cdots,n) \end{array}\right.$$

显然 $\boldsymbol{A}^{(k+1)}$ 的第 $1$ 行直到第 $k$ 行与 $\boldsymbol{A}^{(k)}$ 相同

继续这一过程，直到完成第 $n-1$ 次消元，最后得到与原方程组等价的三角方程组

$$\boldsymbol{A}^{(n)} \boldsymbol{x}=\boldsymbol{b}^{(n)}$$

或

$$\left(\begin{array}{cccc} a_{11}^{(1)} & a_{12}^{(1)} & \cdots & a_{1 n}^{(1)} \\ & a_{22}^{(2)} & \cdots & a_{2 n}^{(2)} \\ & & \ddots & \vdots \\ & & & a_{m n}^{(n)} \end{array}\right)\left(\begin{array}{c} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \end{array}\right)=\left(\begin{array}{c} b_{1}^{(1)} \\ b_{2}^{(2)} \\ \vdots \\ b_{n}^{(n)} \end{array}\right)$$

由原方程组约化为上式的过程称为消元过程。求解三角方程组，设 $a_{i i}^{(i)} \neq 0(i=1,2,\cdots,n-1)$，易得求解公式

$$\left\{\begin{array}{l} x_{n}=b_{n}^{(n)} / a_{n k}^{(n)},\\ x_{k}=\left(b_{k}^{(k)}-\sum_{j=k+1}^{n} a_{k j}^{(k)} x_{j}\right) / a_{k k}^{(k)} \end{array} \quad(k=n-1,n-2,\cdots,2,1)\right.$$

这个求解过程称为回代过程

如果 $a_{11}^{(1)}=0$，那么，由于 $\boldsymbol{A}$ 为非奇异矩阵，所以 $\boldsymbol{A}$ 的第 $1$ 列一定有元素不等于零，例如 $a_{i_{1},1} \neq 0$，于是可交换两行元素 $\left(\boldsymbol{r}_{1} \longleftrightarrow \boldsymbol{r}_{i_{1}}\right)$，将 $a_{i_{1},1}$ 调到第 $1$ 行第 $1$ 列的位置，然后进行消元计算，这时 $\boldsymbol{A}^{(2)}$ 右下角矩阵（$n-1$ 阶）亦为非奇异矩阵。继续这一过程，Gauss 消去法照样可进行计算。

可行性定理

如果 $\boldsymbol{A}$ 为 $n$ 阶非奇异矩阵，则可通过 Gauss 消去法（及交换两行的初等变换）将方程组化为三角方程组。

引理：约化的主元素 $a_{i i}^{(i)} \neq 0(i=1,2,\cdots,k)$ 的充要条件是矩阵 $\boldsymbol{A}$ 的顺序主子式 $D_{i} \neq 0(i=1,2,\cdots,k)$，即

$$\begin{array}{c} D_{1}=a_{11} \neq 0,\\ D_{i}=\left|\begin{array}{ccc} a_{11} & \cdots & a_{1 k} \\ \vdots & & \vdots \\ a_{i 1} & \cdots & a_{i i} \end{array}\right| \neq 0 \quad(i=2,3,\cdots,k) \end{array}$$

计算资源

如果 $\boldsymbol{A}$ 为 $n$ 阶非奇异矩阵，则用 Gauss 消去法解方程组所需乘除法次数及加减法次数分别为

$$\begin{array}{l} \text{M D}=n^{3} / 3+n^{2}-n / 3\\ \text{A S}=n(n-1)(2 n+5) / 6 \end{array}$$

矩阵的三角分解

设系数矩阵 $\boldsymbol{A}$ 的各顺序主子式均不为零，由于对 $\boldsymbol{A}$ 施行行的初等变换相当于用初等矩阵左乘 $\boldsymbol{A}$，于是施行第一次消元后 $\boldsymbol{A}^{(1)}$ 化为 $\boldsymbol{A}^{(2)}$，$\boldsymbol{b}^{(1)}$ 化为 $\boldsymbol{b}^{(2)}$，即

$$\begin{array}{c} \boldsymbol{L}_{1} \boldsymbol{A}^{(1)}=\boldsymbol{A}^{(2)},\quad \boldsymbol{L}_{1} \boldsymbol{b}^{(1)}=\boldsymbol{b}^{(2)},\\ \boldsymbol{L}_{1}=\left(\begin{array}{ccccc} 1 & & & & \\ -m_{21} & 1 & & & \\ -m_{31} & & 1 & & \\ \vdots & & & \ddots & \\ -m_{n 1} & & & & 1 \end{array}\right) \end{array}$$

一般第 $k$ 步消元，$\boldsymbol{A}^{(k)}$ 化为 $\boldsymbol{A}^{(k+1)}$，$\boldsymbol{b}^{(k)}$ 化为 $\boldsymbol{b}^{(k+1)}$，相当于

$$\boldsymbol{L}_{k} \boldsymbol{A}^{(k)}=\boldsymbol{A}^{(k+1)},\quad \boldsymbol{L}_{k} \boldsymbol{b}^{(k)}=\boldsymbol{b}^{(k+1)}$$

重复这一过程，最后得到

$$\left\{\begin{array}{l} \boldsymbol{L}_{n-1} \cdots \boldsymbol{L}_{2} \boldsymbol{L}_{1} \boldsymbol{A}^{(1)}=\boldsymbol{A}^{(n)},\\ \boldsymbol{L}_{n-1} \cdots \boldsymbol{L}_{2} \boldsymbol{L}_{1} \boldsymbol{b}^{(1)}=\boldsymbol{b}^{(n)}, \end{array}\right.$$

其中

$$\boldsymbol{L}_{k}=\left(\begin{array}{cccccc} 1 & & & & & \\ & \ddots & & & & \\ & & 1 & & & \\ & & -m_{k+1,k} & 1 & & \\ & & \vdots & & \ddots & \\ & & -m_{n k} & & & 1 \end{array}\right)$$

将上三角矩阵 $\boldsymbol{A}^{(n)}$ 记作 $\boldsymbol{U}$，于是

$$\boldsymbol{A}=\boldsymbol{L}_{1}^{-1} \boldsymbol{L}_{2}^{-1} \cdots \boldsymbol{L}_{n-1}^{-1} \boldsymbol{U}=\boldsymbol{L} \boldsymbol{U},$$

其中

$$\boldsymbol{L}=\boldsymbol{L}_{1}^{-1} \boldsymbol{L}_{2}^{-1} \cdots \boldsymbol{L}_{n-1}^{-1}=\left(\begin{array}{ccccc} 1 & & & & \\ m_{21} & 1 & & & \\ m_{31} & m_{32} & 1 & & \\ \vdots & \vdots & \ddots & \ddots & \\ m_{n 1} & m_{n 2} & \cdots & m_{n,n-1} & 1 \end{array}\right)$$

为单位下三角矩阵。

这就是说，Gauss 消去法实质上产生了一个将 $\boldsymbol{A}$ 分解为两个三角矩阵相乘的因式分解。

定理：（矩阵的 LU 分解）设 $\boldsymbol{A}$ 为 $n$ 阶矩阵，如果 $\boldsymbol{A}$ 的顺序主子式 $D_{i} \neq 0 ( i =1,2,\cdots,n-1)$，则 $\boldsymbol{A}$ 可分解为一个单位下三角矩阵 $\boldsymbol{L}$ 和一个上三角矩阵 $\boldsymbol{U}$ 的乘积，且这种分解是唯一的。

完全主元素消去法

由 Gauss 消去法知道，在消元过程中可能出现 $a_{k k}^{(k)}=0$ 的情况，这时消去法将无法进行；即使在主元素 $a_{k k}^{(k)} \neq 0$ 但很小时，用其作除数，也会导致其他元素数量级的严重增长和舍入误差的扩散，最后会使得计算解不可靠。

设方程组的增广矩阵为

$$\boldsymbol{B}=\left(\begin{array}{cccccc:c} a_{11} & a_{12} & \cdots & a_{1 j_{1}} & \cdots & a_{1 n} & b_{1} \\ a_{21} & a_{22} & \cdots & a_{2 j_{1}} & \cdots & a_{2 n} & b_{2} \\ \vdots & \vdots & & \vdots & & \vdots & \vdots \\ a_{i_{1} 1} & a_{i_{1} 2} & \cdots & a_{i_{1} j_{1}} & \cdots & a_{i_{1} n} & b_{i_{1}} \\ \vdots & \vdots & & \vdots & & \vdots & \vdots \\ a_{n 1} & a_{n 2} & \cdots & a_{n j_{1}} & \cdots & a_{n n} & b_{n} \end{array}\right)$$

首先在 $\boldsymbol{A}$ 中选取绝对值最大的元素作为主元素，例如 $\left|a_{i_{1} j_{1}}\right|=\max _{\substack{1 \leqslant i \leqslant n \\ 1 \leqslant j \leqslant n}}\left|a_{i j}\right| \neq 0$，然后交换 $\boldsymbol{B}$ 的第 $1$ 行与第 $i_{1}$ 行，第 $1$ 列与第 $j_{1}$ 列，经第一次消元计算得

$$(\boldsymbol{A},\boldsymbol{b}) \rightarrow\left(\boldsymbol{A}^{(2)},\boldsymbol{b}^{(2)}\right)$$

重复上述过程，设已完成第 $k-1$ 步的选主元素，交换两行及交换两列，消元计算，$(\boldsymbol{A},\boldsymbol{b})$ 约化为

$$\left(\boldsymbol{A}^{(k)},\boldsymbol{b}^{(k)}\right)=\left[\begin{array}{cccccc:c} a_{11} & a_{12} & \cdots & \cdots & \cdots & a_{m} & b_{1} \\ & a_{22} & \cdots & \cdots & \cdots & a_{2 n} & b_{2} \\ & & \ddots & & & \vdots & \vdots \\ & & & a_{k k} & \cdots & a_{k n} & b_{k} \\ & & & \vdots & & \vdots & \vdots\\ & & & a_{n k} & \cdots & a_{n n} & b_{n} \end{array}\right]$$

其中 $\boldsymbol{A}^{(k)}$ 元素仍记作 $a_{i j}$，$\boldsymbol{b}^{(k)}$ 元素仍记作 $b_{i}(k=1,2,\cdots,n-1)$
第 $k$ 步选主元素（在 $\boldsymbol{A}^{(k)}$ 右下角方框内选），即确定 $i_{k},j_{k}$ 使

$$\left|a_{i_{k} j_{k}}\right|=\max _{\substack{k \leq i \leq n\\ k \leq j \leq n}}\left|a_{i j}\right| \neq 0$$

交换 $\left(\boldsymbol{A}^{(k)},\boldsymbol{b}^{(k)}\right)$ 第 $k$ 行与 $i_{k}$ 行元素交换 $\left(\boldsymbol{A}^{(k)}\right)$ 第 $k$ 列与 $j_{k}$ 列元素，将 $a_{i_{k} j_{k}}$ 调到 $(k,k)$ 位置，再进行消元计算，最后将原方程组化为

$$\left(\begin{array}{cccc} a_{11} & a_{12} & \cdots & a_{1 n} \\ & a_{22} & \cdots & a_{2 n} \\ & & \ddots & \vdots \\ & & & a_{n n} \end{array}\right)\left(\begin{array}{c} y_{1} \\ y_{2} \\ \vdots \\ y_{n} \end{array}\right)=\left(\begin{array}{c} b_{1} \\ b_{2} \\ \vdots \\ b_{n} \end{array}\right)$$

其中 $y_{1},y_{2},\cdots,y_{n}$ 的次序为未知数 $x_{1},x_{2},\cdots,x_{n}$ 调换后的次序。回代求解得

$$\left\{\begin{array}{l} y_{n}=b_{n} / a_{n n},\\ y_{i}=\left(b_{i}-\sum_{j=i+1}^{n} a_{i j} y_{j}\right) / a_{i i} \quad(i=n-1,\cdots,2,1) \end{array}\right.$$

列主元素消去法

完全主元素消去法在选主元素时要花费较多机器时间。下面介绍另一种常用的方法即列主元素消去法，它仅考虑依次按列选主元素，然后换行使之变到主元位置上，再进行消元计算。

设用列主元素消去法解 $\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}$ 已完成 $k-1$ 步计算，即有

$$(\boldsymbol{A},\boldsymbol{b}) \rightarrow\left(\boldsymbol{A}^{(k)},\boldsymbol{b}^{(k)}\right)=\left(\begin{array}{cccccc：c} a_{11}^{(1)} & a_{12}^{(1)} & \cdots & \cdots & \cdots & a_{1 n}^{(1)} & b_{1}^{(1)} \\ & a_{22}^{(2)} & \cdots & \cdots & \cdots & a_{2 n}^{(2)} & b_{2}^{(2)} \\ & & \ddots & & & \vdots & \vdots \\ & & & a_{k k}^{(k)} & \cdots & a_{k n}^{(k)} & b_{k}^{(k)} \\ & & & \vdots & & \vdots & \vdots \\ & & & a_{n k}^{(k)} & \cdots & a_{n n}^{(k)} & b_{n}^{(n)} \end{array}\right)$$

且 $\boldsymbol{A}^{(k)} x=\boldsymbol{b}^{(k)}$ 与 $\boldsymbol{A x}=\boldsymbol{b}$ 等价，第 $k$ 步选主元素（在 $\boldsymbol{A}^{(k)}$ 第 $k$ 列方框内选），即确定 $i_{k}$ 使

$$\left|a_{i_{k},k}^{(k)}\right|=\max _{k \leqslant i \leqslant n}\left|a_{i k}^{(k)}\right|$$

定理：（列主元素的三角分解定理）如果 $\boldsymbol{A}$ 为非奇异矩阵，则存在排列矩阵 $\boldsymbol{P}$，使

$$\boldsymbol{P} \boldsymbol{A}=\boldsymbol{L} \boldsymbol{U}$$

其中

• $\boldsymbol{L}$ 为单位下三角阵
• $\boldsymbol{U}$ 为上三角阵

直接三角分解法

将 Gauss 消去法改写为紧凑形式，可以直接从矩阵 $\boldsymbol{A}$ 的元素得到计算 $\boldsymbol{L},\boldsymbol{U}$ 元素的递推公式，而不需任何中间步骤，这就是所谓直接三角分解法。一旦实现了矩阵 $\boldsymbol{A}$ 的 $LU$ 分解，那么求解方程组的问题就等价于求解以下两个三角方程组：

• $\boldsymbol{L y}=\boldsymbol{b}$，求 $\boldsymbol{y}$
• $\boldsymbol{U} \boldsymbol{x}=\boldsymbol{y}$，求 $\boldsymbol{x}$

不选主元的三角分解法

设 $\boldsymbol{A}$ 为非奇异矩阵，且有分解式 $\boldsymbol{A}=\boldsymbol{L} \boldsymbol{U}$，其中 $\boldsymbol{L}$ 为单位下三角阵，$\boldsymbol{U}$ 为上三角阵，即

$$\boldsymbol{A}=\left(\begin{array}{cccc} 1 & & & \\ l_{21} & 1 & & \\ \vdots & \ddots & \ddots & \\ l_{n 1} & \cdots & l_{n,n-1} & 1 \end{array}\right)\left(\begin{array}{cccc} u_{11} & u_{12} & \cdots & u_{1 n} \\ & u_{22} & \cdots & u_{2 n} \\ & & \ddots & \vdots \\ & & & u_{n n} \end{array}\right)$$

下面说明 $\boldsymbol{L}$，$\boldsymbol{U}$ 的元素可以由 $n$ 步直接计算定出，其中第 $r$ 步定出 $\boldsymbol{U}$ 的第 $r$ 行和 $\boldsymbol{L}$ 的第 $r$ 列元素。由上式，有

$$a_{1 i}=u_{1 i} \quad(i=1,2,\cdots,n)$$

于是得 $\boldsymbol{U}$ 的第 $1$ 行元素

$$a_{i 1}=l_{i 1} u_{11},\quad l_{i 1}=a_{i 1} / u_{11} \quad(i=2,\cdots,n)$$

于是得 $\boldsymbol{L}$ 的第 $1$ 列元素

设已经定出 $\boldsymbol{U}$ 的第 $1$ 行到第 $r-1$ 行元素与 $\boldsymbol{L}$ 的第 $1$ 列到第 $r-1$ 列元素。利用矩阵乘法，有

$$a_{r i}=\sum_{k=1}^{n} l_{r k} u_{k i}=\sum_{k=1}^{r-1} l_{k k} u_{k i}+u_{r i} \quad\left(\mathrm{当} r<k,l_{r k}=0 \mathrm{时}\right)$$

故

$$u_{r i}=a_{n i}-\sum_{k=1}^{r-1} l_{k_{k}} u_{k i} \quad(i=r,r+1,\cdots,n),$$

又有

$$a_{i r}=\sum_{k=1}^{n} l_{i k} u_{k r}=\sum_{k=1}^{r-1} l_{i k} u_{k r}+l_{i r} u_{r r}$$

总结上述讨论，得到用直接三角分解法解 $\boldsymbol{A x}=b$（要求 $\boldsymbol{A}$ 所有顺序主子式都不为零）的计算公式，步骤如下
步 1：计算 $\boldsymbol{U}$ 的第 $1$ 行元素

$$u_{1 i}=a_{1 i}(i=1,2,\cdots,n),\quad l_{i 1}=a_{i 1} / u_{11}(i=2,3,\cdots,n)$$

步 2：计算 $\boldsymbol{U}$ 第 $r$ 行，$\boldsymbol{L}$ 的第 $r$ 列元素，$r=2,3,\cdots,n$

$$\begin{array}{c} u_{n i}=a_{n i}-\sum_{k=1}^{r-1} l_{r k} u_{k i} \quad(i=r,r+1,\cdots,n) \\ l_{i r}=\left(a_{i r}-\sum_{k=1}^{r-1} l_{i k} u_{k r}\right) / u_{r r} \quad(i=r+1,\cdots,n;r \neq n) \end{array}$$

步 3：求解 $\boldsymbol{L} \boldsymbol{y}=\boldsymbol{b}$

$$\left\{\begin{array}{l} y_{1}=b_{1},\\ y_{i}=b_{i}-\sum_{k=1}^{i-1} l_{i k} y_{k}\quad(i=2,3,\cdots,n) \end{array}\right.$$

步 4：求解 $\boldsymbol{U} \boldsymbol{x}=\boldsymbol{y}$

$$\left\{\begin{array}{l} x_{n}=y_{n} / u_{nn},\\ x_{i}=\left(y_{i}-\sum_{k=i+1}^{n} u_{i k} x_{k}\right) / u_{i i} \quad(i=n-1,n-2,\cdots,1) \end{array}\right.$$

直接分解法大约需要 $n^{3} / 3$ 次乘、除法运算，和 Gauss 消去法的计算量基本相同

选主元的三角分解法

从直接三角分解公式可看出，当 $u_{r r}=0$ 时计算将中断或者当 $u_{r r}$ 绝对值很小时，按分解公式计算可能引起舍入误差的累积。但如果 $\boldsymbol{A}$ 非奇异，就可通过交换 $\boldsymbol{A}$ 的行实现矩阵 $\boldsymbol{P A}$ 的 $LU$ 分解，因此可采用与列主元素消去法类似的方法（可以证明下述方法与列主元素消去法等价），将直接三角分解法修改为（部分）选主元的三角分解法。

（略）

平方根法

应用有限元法解结构力学问题，最后归结为求解线性方程组，这时系数矩阵大多具有对称正定性质。所谓平方根法，就是利用对称正定矩阵的三角分解而得到的求解对称正定方程组的一种有效方法，目前在计算机上广泛应用平方根法解此类方程组。

设 $\boldsymbol{A}$ 为对称阵，且 $\boldsymbol{A}$ 的所有顺序主子式均不为零，则 $\boldsymbol{A}$ 可唯一分解为 $\boldsymbol{A}=\boldsymbol{L} \boldsymbol{U}$ 的形式。为了利用 $\boldsymbol{A}$ 的对称性，将 $\boldsymbol{U}$ 再分解，即

$$\boldsymbol{U}=\left(\begin{array}{llll} u_{11} & & & \\ & u_{22} & & \\ & & \ddots & \\ & & & u_{n n} \end{array}\right)\left(\begin{array}{cccc} 1 & \frac{u_{12}}{u_{11}} & \cdots & \frac{u_{1 n}}{u_{11}} \\ & \ddots & \ddots & \vdots \\ & & \ddots & \frac{u_{n-1,n}}{u_{n-1,n-1}} \\ & & & 1 \end{array}\right)=\boldsymbol{D} \boldsymbol{U}_{0}$$

其中 $\boldsymbol{D}$ 为对角阵，$\boldsymbol{U}_{0}$ 为单位上三角阵 $.$ 于是

$$\boldsymbol{A}=\boldsymbol{L} \boldsymbol{U}=\boldsymbol{L D} \boldsymbol{U}_{0}$$

又

$$\boldsymbol{A}=\boldsymbol{A}^{\mathrm{T}}=\boldsymbol{U}_{0}^{\mathrm{T}}\left(\boldsymbol{D L}^{\mathrm{T}}\right),$$

由分解的唯一性即得 $\boldsymbol{U}_{0}^{\mathrm{T}}=\boldsymbol{L}$，于是得到对称矩阵 $\boldsymbol{A}$ 的分解式 $\boldsymbol{A}=\boldsymbol{L D L}^{\mathrm{T}}$

定理：（对称阵的三角分解定理）设 $\boldsymbol{A}$ 为 $n$ 阶对称阵，且 $\boldsymbol{A}$ 的所有顺序主子式均不为零，则 $\boldsymbol{A}$ 可唯一分解为

$$\boldsymbol{A}=\boldsymbol{L D} \boldsymbol{L}^{\mathrm{T}}$$

其中 $\boldsymbol{L}$ 为单位下三角阵，$\boldsymbol{D}$ 为对角阵。

定理：（对称正定矩阵的三角分解或 Cholesky 分解）如果 $\boldsymbol{A}$ 为 $n$ 阶对称正定矩阵，则存在一个实的非奇异下三角阵 $\boldsymbol{L}$ 使 $\boldsymbol{A}=\boldsymbol{L} \boldsymbol{L}^{\mathrm{T}}$，当限定 $\boldsymbol{L}$ 的对角元素为正时，这种分解是唯一的。

下面用直接分解方法来确定计算 $\boldsymbol{L}$ 元素的递推公式。因为

$$\boldsymbol{A}=\left(\begin{array}{cccc} l_{11} & & & \\ l_{21} & l_{22} & & \\ \vdots & \vdots & \ddots & \\ l_{n 1} & l_{n 2} & \cdots & l_{n n} \end{array}\right)\left(\begin{array}{cccc} l_{11} & l_{21} & \cdots & l_{n 1} \\ & l_{22} & \cdots & l_{n 2} \\ & & \ddots & \vdots \\ & & & l_{n n} \end{array}\right)$$

其中 $l_{i i}>0(i=1,2,\cdots,n)$，由矩阵乘法及 $l_{j k}=0$（当 $j<k$ 时），得

$$a_{i j}=\sum_{k=1}^{n} l_{i k} l_{j k}=\sum_{k=1}^{j-1} l_{i k} l_{j k}+l_{j j} l_{i j}$$

于是得到以下解对称正定方程组 $\boldsymbol{A x}=\boldsymbol{b}$ 的平方根法计算公式

对于 $j=1,2,\cdots,n$，

步 1：计算第 $j$ 个对角线元素

$$l_{j j}=\left(a_{j j}-\sum_{k=1}^{j-1} l_{j k}^{2}\right)^{\frac{1}{2}}$$

步 2：计算下三角矩阵第 $j$ 列所有元素

$$l_{i j}=\left(a_{i j}-\sum_{k=1}^{j-1} l_{i k} l_{j k}\right) / l_{j j} \quad(i=j+1,\cdots,n)$$

步 3：求解 $\boldsymbol{L} \boldsymbol{y}=\boldsymbol{b}$

$$y_{i}=\left(b_{i}-\sum_{k=1}^{i-1} l_{i k} y_{k}\right) / l_{i n} \quad(i=1,2,\cdots,n)$$

步 4：求解 $\boldsymbol{U} \boldsymbol{x}=\boldsymbol{y}$

$$x_{i}=\left(b_{i}-\sum_{k=i+1}^{n} l_{k i} x_{k}\right) / l_{i i} \quad(i=n,n-1,\cdots,1)$$

平方根法约需 $n^{3} / 6$ 次乘除法运算，大约为一般直接 $LU$ 分解法计算量的一半。

追赶法

在一些实际问题中，例如解常微分方程边值问题，解热传导方程以及船体数学放样中建立三次样条函数等中，都会要求解系数矩阵为对角占优的三对角方程组

$$\left(\begin{array}{ccccccc} b_{1} & c_{1} & & & & & \\ a_{2} & b_{2} & c_{2} & & & & \\ & \ddots & \ddots & \ddots & & & \\ & & a_{i} & b_{i} & c_{i} & & \\ & & & \ddots & \ddots & \ddots & \\ & & & & a_{n-1} & b_{n-1} & c_{n-1} \\ & & & & & a_{n} & b_{n} \end{array}\right)\left(\begin{array}{c} x_{1} \\ x_{2} \\ \vdots \\ x_{i} \\ \vdots \\ x_{n-1} \\ x_{n} \end{array}\right)=\left(\begin{array}{c} f_{1} \\ f_{2} \\ \vdots \\ f_{i} \\ \vdots \\ f_{n-1} \\ f_{n} \end{array}\right)$$

简记作

$$\boldsymbol{A x}=\boldsymbol{f}$$

其中 $\boldsymbol{A}$ 满足下列对角占优条件：

• $\left|b_{1}\right|>\left|c_{1}\right|>0$
• $\left|b_{i}\right| \geqslant\left|a_{i}\right|+\left|c_{i}\right|,a_{i},c_{i} \neq 0(i=2,3,\cdots,n-1)$
• $\left|b_{n}\right|>\left|a_{n}\right|>0$

由系数阵 $\boldsymbol{A}$ 的特点，可以将 $\boldsymbol{A}$ 分解为两个三角阵的乘积，即

$$\boldsymbol{A}=\boldsymbol{L} \boldsymbol{U}$$

其中 $\boldsymbol{L}$ 为下三角阵，$\boldsymbol{U}$ 为单位上三角阵。
设

$$\boldsymbol{A}=\left(\begin{array}{ccccc} b_{1} & c_{1} & & & \\ a_{2} & b_{2} & c_{2} & & \\ & \ddots & \ddots & \ddots & \\ & & a_{n-1} & b_{n-1} & c_{n-1} \\ & & & a_{n} & b_{n} \end{array}\right)=\left(\begin{array}{ccccc} \alpha_{1} & & & & \\ \gamma_{2} & \alpha_{2} & & & \\ & \gamma_{3} & \alpha_{3} & & \\ & & \ddots & \ddots & \\ & & & \gamma_{n} & \alpha_{n} \end{array}\right)\left(\begin{array}{ccccc} 1 & \beta_{1} & & & \\ & 1 & \beta_{2} & & \\ & & \ddots & \ddots & \\ & & & 1 & \beta_{n-1} \\ & & & & 1 \end{array}\right)$$

其中 $\alpha_{i},\beta_{i},\gamma_{i}$ 为待定系数，比较上式两端即得

$$\left.\begin{array}{ll} b_{1}=\alpha_{1},& c_{1}=\alpha_{1} \beta_{1},\\ a_{i}=\gamma_{i},& b_{i}=\gamma_{i} \beta_{i-1}+\alpha_{i} \quad(i=2,\cdots,n),\\ c_{i}=\alpha_{i} \beta_{i} & (i=2,3,\cdots,n-1) \end{array}\right\}$$

步 1：计算 $\left\{\beta_{i}\right\}$ 的递推公式

$$\beta_{1}=c_{1} / b_{1},\quad \beta_{i}=c_{i} /\left(b_{i}-a_{i} \beta_{i-1}\right) \quad(i=2,3,\cdots,n-1)$$

步 2：解 $\boldsymbol{L} \boldsymbol{y}=\boldsymbol{f}$

$$\left\{\begin{array}{l} y_{1}=f_{1} / b_{1}\\ y_{i}=\left(f_{i}-a_{i} y_{i-1}\right) /\left(b_{i}-a_{i} \beta_{i-1}\right) \quad(i=2,3,\cdots,n) \end{array}\right.$$

步 3：解 $\boldsymbol{U} \boldsymbol{x}=\boldsymbol{y}$

$$\left\{\begin{array}{l} x_{n}=y_{n}\\ x_{i}=y_{i}-\beta_{i} x_{i+1} \quad(i=n-1,n-2,\cdots,2,1) \end{array}\right.$$

范数

数量积

我们用 $\mathbf{R}^{n}$ 表示 $n$ 维实向量空间，用 $\mathbf{C}^{n}$ 表示 $n$ 维复向量空间，首先将向量长度概念推广到 $\mathbf{R}^{n}$（或 $\mathbf{C}^{n}$）中：
设

$$\boldsymbol{x}=\left(x_{1},x_{2},\cdots,x_{n}\right)^{\mathrm{T}},\quad \boldsymbol{y}=\left(y_{1},y_{2},\cdots,y_{n}\right)^{\mathrm{T}} \in \mathbf{R}^{n} \left(\mathrm{或} \mathbf{C}^{n}\right)$$

将实数 $(\boldsymbol{x},\boldsymbol{y})=\boldsymbol{y}^{\mathrm{T}} \boldsymbol{x}=\sum_{i=1}^{n} x_{i} y_{i}$ （或复数 $(\boldsymbol{x},\boldsymbol{y})=\boldsymbol{y}^{H} \boldsymbol{x}=\sum_{i=1}^{n} x_{i} \overline{y_{i}}$，其中 $\boldsymbol{y}^{H}=\overline{\boldsymbol{y}}^{\mathrm{T}}$ ）称为向量 $\boldsymbol{x},\boldsymbol{y}$ 的数量积

欧拉范数

非负实数

$$\|\boldsymbol{x}\|_{2}=(\boldsymbol{x},\boldsymbol{x})^{\frac{1}{2}}=\left(\sum_{i=1}^{n} x_{i}^{2}\right)^{\frac{1}{2}}$$

称为向量 $\boldsymbol{x}$ 的 Euclid 范数

数量积的性质

设 $x,y \in \mathbf{R}^{n}$ （或 $\mathbf{C}^{n}$），则

• $(\boldsymbol{x},\boldsymbol{x})=0$，当且仅当 $\boldsymbol{x}=\mathbf{0}$ 时成立
• $(\alpha \boldsymbol{x},\boldsymbol{y})=\alpha(\boldsymbol{x},\boldsymbol{y}),\alpha$ 为实数（或 $(\boldsymbol{x},\alpha \boldsymbol{y})=\bar{\alpha}(\boldsymbol{x},\boldsymbol{y}),\alpha$ 为复数）
• $(\boldsymbol{x},\boldsymbol{y})=(\boldsymbol{y},\boldsymbol{x})$ （或 $(\boldsymbol{x},\boldsymbol{y})=(\overline{\boldsymbol{y},\boldsymbol{x}})$
• $\left(\boldsymbol{x}_{1}+\boldsymbol{x}_{2},\boldsymbol{y}\right)=\left(\boldsymbol{x}_{1},\boldsymbol{y}\right)+\left(\boldsymbol{x}_{2},\boldsymbol{y}\right)$
• （Cauchy-Schwarz 不等式）$|(\boldsymbol{x},\boldsymbol{y})| \leqslant|\boldsymbol{x}|_{2} \cdot\|\boldsymbol{y}\|_{2}$，等式当且仅当 $x$ 与 $y$ 线性相关时成立
• （三角不等式）$\|\boldsymbol{x}+\boldsymbol{y}\|_{2} \leqslant\|\boldsymbol{x}\|_{2}+\|\boldsymbol{y}\|_{2}$

向量范数的一般定义

如果向量 $\boldsymbol{x} \in \mathbf{R}^{n}$（或 $\mathbf{C}^{n}$）的某个实值函数 $N(\boldsymbol{x})=\|x\|$，满足条件

1. $\|\boldsymbol{x}\| \geqslant 0$（$\|\boldsymbol{x}\|=0$ 当且仅当 $\boldsymbol{x}=\mathbf{0}$）（正定条件）
2. $\|\alpha \boldsymbol{x}\|=|\alpha|\|\boldsymbol{x}\|,\forall \alpha \in \mathbf{R}$（或 $\alpha \in \mathbf{C}$）
3. $\|\boldsymbol{x}+\boldsymbol{y}\| \leqslant\|\boldsymbol{x}\|+\|\boldsymbol{y}\|$（三角不等式）

则称 $N(\boldsymbol{x})$ 是 $\mathbf{R}^{n}$（或 $\mathbf{C}^{n}$）上的一个向量范数（或模）

• 由条件 3 可推出不等式 $|\|\boldsymbol{x}\|-\|\boldsymbol{y}\|| \leqslant\|\boldsymbol{x}-\boldsymbol{y}\|$

常用的向量范数

• 向量的 $\infty -$ 范数（最大范数）： $$\|\boldsymbol{x}\|_{\infty}=\max _{1 \leqslant i \leqslant n}\left|x_{i}\right|$$
• 向量的 $1-$ 范数： $$\|\boldsymbol{x}\|_{1}=\sum_{i=1}^{n}\left|x_{i}\right|$$
• 向量的 $2-$ 范数： $$\|\boldsymbol{x}\|_{2}=(\boldsymbol{x},\boldsymbol{x})^{\frac{1}{2}}=\left(\sum_{i=1}^{n} x_{i}^{2}\right)^{\frac{1}{2}}$$ $N(\boldsymbol{x})= |\boldsymbol{x}|_{2}$ 是 $\mathbf{R}^{n}$ 上一个向量范数，称为向量 $\boldsymbol{x}$ 的 Euclid 范数
• 向量的 $p-$ 范数： $$\|\boldsymbol{x}\|_{p}=\left(\sum_{i=1}^{n}\left|x_{i}\right|^{p}\right)^{1 / p},\quad p \in[1,\infty)$$

向量范数的性质

• （$N(\boldsymbol{x})$ 的连续性）设非负函数 $N(\boldsymbol{x})=\|\boldsymbol{x}\|$ 为 $\mathbf{R}^{n}$ 上任一向量范数，则 $N(\boldsymbol{x})$ 是 $\boldsymbol{x}$ 分量 $x_{1},x_{2},\cdots,x_{n}$ 的连续函数。
• （向量范数的等价性）设 $\|\boldsymbol{x}\|_{s},\|\boldsymbol{x}\|_{t}$，为 $\mathbf{R}^{n}$ 上向量的任意两种范数，则存在常数 $c_{1},c_{2}>0$，使得 $$c_{1}\|\boldsymbol{x}\|_{s} \leqslant\|\boldsymbol{x}\|_{t} \leqslant c_{2}\|\boldsymbol{x}\|_{s} \mathrm{，对一切} \boldsymbol{x} \in \mathbf{R}^{n}$$
• 如果在一种范数意义下向量序列收敛，则在任何一种范数意义下该向量序列亦收敛
• $\lim _{k \rightarrow \infty} x^{(k)}=x^{*} \Longleftrightarrow \|\boldsymbol{x}^{(k)}-\boldsymbol{x}^{*} \| \rightarrow 0$ （当 $k \rightarrow \infty$ 时），其中 $\|\cdot\|$ 为向量的任一种范数

矩阵范数

如果矩阵 $\boldsymbol{A} \in \mathbf{R}^{n \times n}$ 的某个非负的实值函数 $N(\boldsymbol{A})= \|\boldsymbol{A}\|$，满足条件

1. 正定条件：$\|\boldsymbol{A}\| \geqslant 0(\|\boldsymbol{A}\|=0 \Leftrightarrow \boldsymbol{A}=0)$
2. 齐次条件：$\|c \boldsymbol{A}\|=|c|\|\boldsymbol{A}\|$，$c$ 为实数
3. 三角不等式：$\|\boldsymbol{A}+\boldsymbol{B}\| \leqslant\|\boldsymbol{A}\|+\|\boldsymbol{B}\|$
4. $\|\boldsymbol{A} \boldsymbol{B}\| \leqslant\|\boldsymbol{A}\|\|\boldsymbol{B}\|$

则称 $N(\boldsymbol{A})$ 是 $\mathbf{R}^{n \times n}$ 上的一个矩阵范数（或模）

矩阵的算子范数

设 $\boldsymbol{x} \in \mathbf{R}^{n}$，$\boldsymbol{A} \in \mathbf{R}^{n \times n}$，给出一种向量范数 $\|\boldsymbol{x}\|_{v}$（如 $v=1,2$ 或 $\infty$），相应地定义一个矩阵的算子范数为

$$\|\boldsymbol{A}\|_{v}=\max _{x \neq 0} \frac{\|\boldsymbol{A} \boldsymbol{x}\|_{v}}{\|\boldsymbol{x}\|_{v}}$$

设 $\|\boldsymbol{x}\|_{v}$ 是 $\mathbf{R}^{n}$ 上一个向量范数，则 $\|\boldsymbol{A}\|_{v}$ 是 $\mathbf{R}^{n \times n}$ 上的矩阵范数，且满足相容条件

$$\|\boldsymbol{A} \boldsymbol{x}\|_{v} \leqslant\|\boldsymbol{A}\|_{v}\|\boldsymbol{x}\|_{v}$$

常用矩阵范数

设 $x \in \mathbf{R}^{n}$，$A \in \mathbf{R}^{n \times n}$，则

• $\boldsymbol{A}$ 的行范数：$\|\boldsymbol{A}\|_{\infty}=\max _{1 \leqslant i \leqslant n} \sum_{j=1}^{n}\left|a_{i j}\right|$
• $\boldsymbol{A}$ 的列范数：$\|\boldsymbol{A}\|_{1}=\max _{1 \leqslant j \leqslant n} \sum_{i=1}^{n}\left|a_{i j}\right|$
• $\boldsymbol{A}$ 的 $2-$ 范数：$\|\boldsymbol{A}\|_{2}=\sqrt{\lambda_{\mathrm{max}}\left(\boldsymbol{A}^{\mathrm{T}} \boldsymbol{A}\right)}$

谱半径

设 $\boldsymbol{A} \in \mathbf{R}^{n \times n}$ 的特征值为 $\lambda_{i}(i=1,2,\cdots,n)$，称 $\rho(\boldsymbol{A})=\max _{1 \leqslant i \leqslant n}\left|\lambda_{i}\right|$ 为 $\boldsymbol{A}$ 的谱半径

（特征值上界）设 $\boldsymbol{A} \in \mathrm{R}^{n \times n}$，则 $\rho(\boldsymbol{A}) \leqslant\|\boldsymbol{A}\|$，即 $\boldsymbol{A}$ 的谱半径不超过 $\boldsymbol{A}$ 的任何一种算子范数（对于 $\|\boldsymbol{A}\|_{F}$ 亦对）。

其他内容（略）

如果 $\boldsymbol{A} \in \mathbf{R}^{n \times n}$ 为对称矩阵，则 $\|\boldsymbol{A}\|_{2}=\rho(\boldsymbol{A})$

如果 $\|\boldsymbol{B}\|<1$，则 $\boldsymbol{I} \pm \boldsymbol{B}$ 为非奇异矩阵，且 $\left\|(\boldsymbol{I} \pm \boldsymbol{B})^{-1}\right\| \leqslant \frac{1}{1-\|\boldsymbol{B}\|}$，其中 $\|\cdot\|$ 是指矩阵的算子范数。

病态矩阵

如果矩阵 $\boldsymbol{A}$ 或常数项 $\boldsymbol{b}$ 的微小变化，引起方程组 $\boldsymbol{A x}=\boldsymbol{b}$ 解的巨大变化，则称此方程组为病态方程组，矩阵 $\boldsymbol{A}$ 称为病态矩阵（相对于方程组而言），否则称方程组为良态方程组，$\boldsymbol{A}$ 称为良态矩阵。

相对误差上界

设 $\boldsymbol{A}$ 是非奇异阵，$\boldsymbol{A x}=\boldsymbol{b} \neq \mathbf{0}$，且 $\boldsymbol{A}(\boldsymbol{x}+\delta \boldsymbol{x})=\boldsymbol{b}+\delta \boldsymbol{b}$，则

$$\frac{\|\delta \boldsymbol{x}\|}{\|\boldsymbol{x}\|} \leqslant\left\|\boldsymbol{A}^{-1}\right\|\|\boldsymbol{A}\| \frac{\|\delta \boldsymbol{b}\|}{\|\boldsymbol{b}\|}$$

上式给出了解的相对误差的上界，常数项 $\boldsymbol{b}$ 的相对误差在解中可能放大 $\left\|\boldsymbol{A}^{-1}\right\|\|\boldsymbol{A}\|$ 倍。

设 $\boldsymbol{A}$ 为非奇异矩阵，$\boldsymbol{A} \boldsymbol{x}=\boldsymbol{b} \neq \mathbf{0}$，且 $(\boldsymbol{A}+\delta \boldsymbol{A})(\boldsymbol{x}+\delta \boldsymbol{x})=\boldsymbol{b}$，如果 $\left\|\boldsymbol{A}^{-1}\right\|\|\delta \boldsymbol{A}\|<1$，则

$$\frac{\|\delta \boldsymbol{x}\|}{\|\boldsymbol{x}\|} \leqslant \frac{\left\|\boldsymbol{A}^{-1}\right\|\|\boldsymbol{A}\| \frac{\|\delta \boldsymbol{A}\|}{\|\boldsymbol{A}\|}}{1-\left\|\boldsymbol{A}^{-1}\right\|\left\|\boldsymbol{A}\right\| \frac{\|\delta \boldsymbol{A}\|}{\|\boldsymbol{A}\|}}$$

如果 $\delta \boldsymbol{A}$ 充分小，且在 $\left\|\boldsymbol{A}^{-1}\right\|\|\delta \boldsymbol{A}\|<1$ 条件下，那么上式说明矩阵 $\boldsymbol{A}$ 的相对误差 $\frac{\|\delta \boldsymbol{A}\|}{\|\boldsymbol{A}\|}$ 在解中可能放大 $\left\|\boldsymbol{A}^{-1}\right\|\|\boldsymbol{A}\|$ 倍。

条件数

设 $\boldsymbol{A}$ 为非奇异矩阵，称数 $\operatorname{cond}(\boldsymbol{A})_{v}=\left\|\boldsymbol{A}^{-1}\right\|_{v}\|\boldsymbol{A}\|_{v}$（$v=1,2$ 或 $\infty$）为矩阵 $\boldsymbol{A}$ 的条件数。

• 当 $\boldsymbol{A}$ 的条件数相对地大时，即 $\operatorname{cond}(\boldsymbol{A}) \gg 1$ 时，方程组是病态的（即 $\boldsymbol{A}$ 是病态矩阵，或者说 $\boldsymbol{A}$ 是坏条件的）
• 当 $\boldsymbol{A}$ 的条件数相对地小时，方程组是良态的（或者说 $\boldsymbol{A}$ 是好条件的）
• $\boldsymbol{A}$ 的条件数愈大，方程组的病态程度愈严重，也就愈难得到方程组的比较准确的解

通常使用的条件数有

1. $$\operatorname{cond}(\boldsymbol{A})_{\infty}=\left\|\boldsymbol{A}^{-1}\right\|_{\infty}\|\boldsymbol{A}\|_{\infty}$$
2. $\boldsymbol{A}$ 的谱条件数 $$\operatorname{cond}(\boldsymbol{A})_{2}=\|\boldsymbol{A}\|_{2}\left\|\boldsymbol{A}^{-1}\right\|_{2}=\sqrt{\frac{\lambda_{\max }\left(\boldsymbol{A}^{\mathrm{T}} \boldsymbol{A}\right)}{\lambda_{\min }\left(\boldsymbol{A}^{\mathrm{T}} \boldsymbol{A}\right)}}$$ 当 $\boldsymbol{A}$ 为对称阵时，$\operatorname{cond}(\boldsymbol{A})_{2}=\frac{\left|\lambda_{1}\right|}{\left|\lambda_{n}\right|}$，其中 $\lambda_{1},\lambda_{n}$ 为 $\boldsymbol{A}$ 的绝对值最大和绝对值最小的特征值。

性质

• 对任何非奇异矩阵 $\boldsymbol{A}$，都有 $\operatorname{cond}(\boldsymbol{A})_{v} \geqslant 1$，事实上， $$\operatorname{cond}(\boldsymbol{A})_{v}=\left\|\boldsymbol{A}^{-1}\right\|_{v}\|\boldsymbol{A}\|_{v} \geqslant\left\|\boldsymbol{A}^{-1} \boldsymbol{A}\right\|_{v}=1$$
• 设 $\boldsymbol{A}$ 为非奇异矩阵，$c$ 为不等于零的常数，则 $\operatorname{cond}(c \boldsymbol{A})_{v}=\operatorname{cond}(\boldsymbol{A})_{v}$。
• 如果 $\boldsymbol{A}$ 为正交矩阵，则 $\operatorname{cond}(\boldsymbol{A})_{2}=1$；如果 $\boldsymbol{A}$ 为非奇异矩阵，$\boldsymbol{R}$ 为正交矩阵，
  则 $$\operatorname{cond}(\boldsymbol{R} \boldsymbol{A})_{2}=\operatorname{cond}(\boldsymbol{A R})_{2}=\operatorname{cond}(\boldsymbol{A})_{2}$$

第八章 解线性方程组的迭代法

背景

对于有很多零元素的句子，在计算机内存和运算两方面，迭代法通常都可利用 $\boldsymbol{A}$ 中有大量零元素的特点。

定义

对于给定方程组 $\boldsymbol{x}=\boldsymbol{B} \boldsymbol{x}+\boldsymbol{f}$，设有唯一解 $\boldsymbol{x}^{*}$，则

$$\boldsymbol{x}^{*}=\boldsymbol{B} \boldsymbol{x}^{*}+\boldsymbol{f}$$

又设 $\boldsymbol{x}^{(0)}$ 为任取的初始向量，按下述公式构造向量序列

$$\left\{\begin{array}{c} \boldsymbol{x}^{(1)}=\boldsymbol{B} \boldsymbol{x}^{(0)}+\boldsymbol{f},\\ \boldsymbol{x}^{(2)}=\boldsymbol{B} x^{(1)}+\boldsymbol{f},\\ \vdots \\ \boldsymbol{x}^{(k+1)}=\boldsymbol{B} x^{(k)}+\boldsymbol{f}, \end{array}\right.$$

其中 $k$ 表示迭代次数。

对于给定的方程组 $\boldsymbol{x}=\boldsymbol{B} \boldsymbol{x}+\boldsymbol{f}$，用上式逐步代入求近似解的方法称为迭代法（或称为一阶定常迭代法，这里 $\boldsymbol{B}$ 与 $k$ 无关）

如果 $\lim _{k \rightarrow \infty} \boldsymbol{x}^{(k)}$ 存在（记作 $\boldsymbol{x}^{*}$），称此迭代法收敛，显然 $\boldsymbol{x}^{*}$ 就是方程组的解，否则称此迭代法发散。

引进误差向量

$$\boldsymbol{\varepsilon}^{(k+1)}=\boldsymbol{x}^{(k+1)}-\boldsymbol{x}^{*}$$

由条件得

$$\boldsymbol{\varepsilon}^{(k+1)}=\boldsymbol{B}^{(k)} \quad(k=0,1,2 \cdots)$$

递推得

$$\boldsymbol{\varepsilon}^{(k)}=\boldsymbol{B}^{(k-1)}=\cdots=\boldsymbol{B}^{k} \boldsymbol{\varepsilon}^{(0)}$$

要考察 $\left\{\boldsymbol{x}^{(k)}\right\}$ 的收敛性，就要研究 $\boldsymbol{B}$ 在什么条件下有 $\boldsymbol{\varepsilon}^{(k)} \rightarrow 0(k \rightarrow \infty)$，亦即要研究 $\boldsymbol{B}$ 满足什么条件时有 $\lim _{k \rightarrow \infty} \boldsymbol{B}^{k} \rightarrow \boldsymbol{O}$（零矩阵）

Jacobi 迭代法

设有方程组

$$\sum_{j=1}^{n} a_{i j} x_{j}=b_{i} \quad(i=1,2,\cdots,n)$$

记作

$$\boldsymbol{A x}=\boldsymbol{b}$$

$\boldsymbol{A}$ 为非奇异阵且 $a_{i j} \neq 0(i=1,2,\cdots,n)$。将 $\boldsymbol{A}$ 分裂为 $\boldsymbol{A}=\boldsymbol{D}-\boldsymbol{L}-\boldsymbol{U}$，其中

$$\begin{array}{c} \boldsymbol{D}=\left(\begin{array}{ccccc} a_{11} & & & & \\ & a_{22} & & & \\ & & \ddots & & \\ & & & \ddots & \\ & & & & a_{n n} \end{array}\right)\\ \boldsymbol{L}=-\left(\begin{array}{ccccc} 0 & & & \\ a_{21} & 0 & & \\ a_{31} & a_{32} & 0 & \\ \vdots & \vdots & \ddots & \ddots & \\ a_{n 1} & a_{n 2} & \cdots & a_{n,n-1} & 0 \end{array}\right)\\ \boldsymbol{U}=-\left(\begin{array}{ccccc} 0 & a_{12} & a_{13} & \cdots & a_{1 n} \\ & 0 & a_{23} & \cdots & a_{2 n} \\ & & \ddots & \ddots & \vdots \\ & & & 0 & a_{n-1,n} \\ & & & & 0 \end{array}\right) \end{array}$$

移项整理，得到等价方程组

$$x_{i}=\frac{1}{a_{i i}}\left(b_{i}-\sum_{\substack{j=1 \\ j \neq i}}^{n} a_{i j} x_{j}\right) \quad(i=1,2,\cdots,n),$$

简记作

$$\boldsymbol{x}=\boldsymbol{B}_{0} \boldsymbol{x}+\boldsymbol{f}$$

其中

$$\boldsymbol{B}_{0}=\boldsymbol{I}-\boldsymbol{D}^{-1} \boldsymbol{A}=\boldsymbol{D}^{-1}(\boldsymbol{L}+\boldsymbol{U}),\quad \boldsymbol{f}=\boldsymbol{D}^{-1} \boldsymbol{b}$$

Jacobi 迭代公式

$$\left\{\begin{array}{l} \boldsymbol{x}^{(0)}=\left(x_{1}^{(0)},x_{2}^{(0)},\cdots,x_{n}^{(0)}\right)^{\mathrm{T}} \quad \mathrm{(初始向量)}\\ x_{i}^{(k+1)}=\frac{1}{a_{i i}}\left(b_{i}-\sum_{\substack{i=1 \\ j \neq i}}^{n} a_{i j} x_{j}^{(k)}\right) \end{array}\right.$$

• 其中 $\boldsymbol{x}^{(k)}=\left(x_{1}^{(k)},x_{2}^{(k)},\cdots,x_{n}^{(k)}\right)^{\mathrm{T}}$ 为第 $k$ 次迭代向量。
• 设 $\boldsymbol{x}^{(k)}$ 已经算出，可计算下一次迭代向量 $\boldsymbol{x}^{(k+1)}(k=0,1,2,\cdots;i=1,2,\cdots,n)$
• 迭代公式阵形式为 $$\left\{\begin{array}{l} \boldsymbol{x}^{(0)} \quad \mathrm{(初始向量)}\\ \boldsymbol{x}^{(k+1)}=\boldsymbol{B}_{0} \boldsymbol{x}^{(k)}+\boldsymbol{f} \end{array}\right.$$ 其中 $\boldsymbol{B}_{0}$ 称为 Jacobi 方法迭代矩阵

由此看出，Jacobi 迭代法公式简单，每迭代一次只需计算一次矩阵和向量乘法。在用计算机计算时，需要两组工作单元，以存储 $\boldsymbol{x}^{(k)}$ 及 $\boldsymbol{x}^{(k+1)}$。

Gauss-Seidel 迭代法

对这些最新计算出来的第 $k+1$ 次近似 $\boldsymbol{x}^{(k+1)}$ 的分量 $x_{j}^{(k+1)}$ 加以利用，就得到所谓解方程组的 Gauss-Seidel 迭代法（简称 G-S 方法）

$$\begin{array}{c} x^{(0)}=\left(x_{1}^{(0)},x_{2}^{(0)},\cdots,x_{n}^{(0)}\right)^{\mathrm{T}} \quad\mathrm{(初始向量)}\\ x_{i}^{(k+1)}=\frac{1}{a_{i i}}\left(b_{i}-\sum_{j=1}^{i-1} a_{i j} x_{j}^{(k+1)}-\sum_{j=i+1}^{n} a_{i j} x_{j}^{(k)}\right) \quad(k=0,1,2,\cdots;i=1,2,\cdots,n) \end{array}$$

或写为

$$\left\{\begin{array}{l} x_{i}^{(k+1)}=x_{i}^{(k)}+\Delta x_{i} \quad(k=0,1,2,\cdots;i=1,2,\cdots,n)\\ \Delta x_{i}=\frac{1}{a_{i i}}\left(b_{i}-\sum_{j=1}^{i-1} a_{i j} x_{i}^{(k+1)}-\sum_{j=i}^{n} a_{i j} x_{j}^{(k)}\right) \end{array}\right.$$

矩阵形式

$$\begin{array}{c} \boldsymbol{x}^{(k+1)}=\boldsymbol{G} \boldsymbol{x}^{(k)}+\boldsymbol{f}\\ \boldsymbol{G}=(\boldsymbol{D}-\boldsymbol{L})^{-1} \boldsymbol{U},\quad \boldsymbol{f}=(\boldsymbol{D}-\boldsymbol{L})^{-1} \boldsymbol{b} \end{array}$$

• Gauss-Seidel 迭代法比 Jacobi 迭代法收敛快
• 但有的方程组，Jacobi 方法收敛，而 Gauss-Seidel 迭代法却是发散的

收敛性判断

矩阵序列收敛

设有矩阵序列 $\boldsymbol{A}_{k}=\left(a_{i j}^{(k)}\right)_{n \times n}(k=1,2,\cdots)$ 及 $\boldsymbol{A}=\left(a_{i j}\right)_{n \times n}$，如果

$$\lim _{k \rightarrow \infty} a_{i j}^{(k)}=a_{i j} \quad(i,j=1,2,\cdots,n)$$

成立，则称 $\left\{\boldsymbol{A}_{k}\right\}$ 收敛于 $\boldsymbol{A}$，记作 $\lim _{k \rightarrow \infty} \boldsymbol{A}_{k}=\boldsymbol{A}$。

$\lim_{k \rightarrow \infty} \boldsymbol{A}_{k}=\boldsymbol{A}$ 的充要条件是 $\left|\boldsymbol{A}_{k}-\boldsymbol{A}\right| \rightarrow 0(k \rightarrow \infty)$。

迭代法收敛性定理

设 $\boldsymbol{B}=\left(b_{i j}\right)_{n \times n}$，则 $\boldsymbol{B}^{k} \rightarrow \boldsymbol{O}(k \rightarrow \infty)$ 的充要条件是 $\rho(\boldsymbol{B})<1$。

（迭代法基本定理）设有方程组

$$\boldsymbol{x}=\boldsymbol{B} \boldsymbol{x}+\boldsymbol{f}$$

对于任意初始向量 $\boldsymbol{x}^{(0)}$ 及任意 $\boldsymbol{f}$，解此方程组的迭代法（即 $\boldsymbol{x}^{(k+1)}=\boldsymbol{B} \boldsymbol{x}^{(k)}+\boldsymbol{f}$）收敛的充要条件是 $\rho(\boldsymbol{B})<1$。

（迭代法收敛的充分条件）如果方程组的迭代公式为 $\boldsymbol{x}^{(k+1)}= \boldsymbol{B} \boldsymbol{x}^{(k)}+\boldsymbol{f}$（$\boldsymbol{x}^{(0)}$ 为任意初始向量），且迭代矩阵的某一种范数 $\|\boldsymbol{B}\|_{v}=q<1$，则：

1. 迭代法收敛
2. $$\left\|\boldsymbol{x}^{*}-\boldsymbol{x}^{(k)}\right\|_{v} \leqslant \frac{q}{1-q}\left\|\boldsymbol{x}^{(k)}-\boldsymbol{x}^{(k-1)}\right\|_{v}$$
3. $$\left\|\boldsymbol{x}^{*}-\boldsymbol{x}^{(k)}\right\|_{v} \leqslant \frac{q^{k}}{1-q}\left\|\boldsymbol{x}^{(1)}-\boldsymbol{x}^{(0)}\right\|_{v}$$

Gauss-Seidel 迭代法收敛的充要条件是 $\rho(\boldsymbol{G})<1$，其中 $\boldsymbol{G}$ 为 Gauss-Seidel 迭代法的迭代矩阵。

迭代法的速度

考察误差向量 $\boldsymbol{\varepsilon}^{(k)}=\boldsymbol{x}^{(k)}-\boldsymbol{x}^{*}=\boldsymbol{B}^{k} \boldsymbol{\varepsilon}^{(0)}$。设 $\boldsymbol{B}$ 有 $n$ 个线性无关的特征向量 $\boldsymbol{u}_{1},\boldsymbol{u}_{2},\cdots,\boldsymbol{u}_{n}$，相应的特征值为 $\lambda_{1},\lambda_{2},\cdots,\lambda_{n}$。由 $\boldsymbol{\varepsilon}^{(0)}=\sum_{i=1}^{n} a_{i} \boldsymbol{u}_{i}$，得

$$\boldsymbol{\varepsilon}^{(k)}=\boldsymbol{B}^{k} \boldsymbol{\varepsilon}^{(0)}=\sum_{i=1}^{n} a_{i} \boldsymbol{B}^{k} \boldsymbol{u}_{i}=\sum_{i=1}^{n} a_{i} \lambda_{i}^{k} \boldsymbol{u}_{i}$$

可以看出，当 $\rho(\boldsymbol{B})<1$ 愈小时，$\lambda_{i}^{k} \rightarrow 0(i=1,2,\cdots,n;k \rightarrow \infty)$ 愈快，即 $\boldsymbol{\varepsilon}^{(k)} \rightarrow \mathbf{0}$ 愈快，故可用量 $\rho(\boldsymbol{B})$ 来刻画迭代法的收敛快慢。现在依据给定精度要求来确定迭代次数 $k$，即使

$$[\rho(\boldsymbol{B})]^{k} \leqslant 10^{-s}$$

取对数得

$$k \geqslant \frac{s \ln 10}{-\ln \rho(\boldsymbol{B})}$$

称 $R(B)=-\ln \rho(B)$ 为迭代法的收敛速度。

由于一般当 $n$ 较大时，矩阵特征值计算比较困难，基本定理的条件比较难验证，所以最好建立与矩阵元素直接有关的条件来判别迭代法的收敛性。由于 $\rho(\boldsymbol{B}) \leqslant \|\boldsymbol{B}\|_{v}(v=1,2,\infty \text{或} F)$，所以可用 $\|\boldsymbol{B}\|_{v}$ 来作为 $\rho(\boldsymbol{B})$ 上界的一种估计。

特殊方程的收敛

对角占优阵

设 $\boldsymbol{A}=\left(a_{i j}\right)_{n \times n} \in \mathbf{R}^{n \times n}$（或 $\mathbf{C}^{n \times n}$）

• 如果矩阵 $\boldsymbol{A}$ 满足条件 $$\left|a_{i i}\right|>\sum_{\substack{j=1 \\ j \neq i}}^{n}\left|a_{i j}\right| \quad(i=1,2,\cdots,n)$$ 即 $\boldsymbol{A}$ 的每一行对角元素的绝对值都严格大于同行其他元素绝对值之和，则称 $\boldsymbol{A}$ 为严格对角占优阵。
• 如果 $\left|a_{i i}\right| \geqslant \sum_{\substack{j=1 \ j \neq i}}^{n}\left|a_{i j}\right|(i=1,2,\cdots,n)$ 且至少有一个不等式严格成立，称 $\boldsymbol{A}$ 为弱对角占优阵

（对角占优定理）如果 $\boldsymbol{A}=\left(a_{i j}\right)_{n} \in \mathbf{R}^{n \times n}$（或 $\mathbf{C}^{n \times n}$）为严格对角占优阵或为不可约弱对角占优阵，则 $\boldsymbol{A}$ 是非奇异矩阵

如果 $\boldsymbol{A} \in \mathbf{R}^{n \times m}$ 为严格对角占优阵或为不可约弱对角占优阵，则对于任意的 $\boldsymbol{x}^{(0)}$，解方程组的 Jacobi 迭代法，Gauss-Seidel 迭代法均收敛。

可约与不可约矩阵

设 $\boldsymbol{A}=\left(a_{i j}\right)_{n} \in \mathbf{R}^{n \times n}$（或 $\mathbf{C}^{n \times n}$），当 $n \geqslant 2$ 时，如果存在 $n$ 阶置换阵 $\boldsymbol{P}$ 使

$$\boldsymbol{P}^{\mathrm{T}} \boldsymbol{A} \boldsymbol{P}=\left(\begin{array}{cc} \boldsymbol{A}_{11} & \boldsymbol{A}_{12} \\ \boldsymbol{O} & \boldsymbol{A}_{22} \end{array}\right)$$

成立，其中 $\boldsymbol{A}_{11}$ 为 $r$ 阶子矩阵，$\boldsymbol{A}_{22}$ 为 $n-r$ 阶子矩阵 $(1 \leqslant r \leqslant n)$，则称 $\boldsymbol{A}$ 是可约矩阵。如果不存在置换阵 $\boldsymbol{P}$ 使上式成立，则称 $\boldsymbol{A}$ 是不可约矩阵。

$\boldsymbol{A}$ 是可约矩阵，意味着 $\boldsymbol{A x}=\boldsymbol{b}$ 可经过若干行列重排（若 $\boldsymbol{A}$ 经过两行交换的同时进行相应的两列的交换，称对 $\boldsymbol{A}$ 进行一次行列重排），化为两个低阶方程组求解。

事实上，由 $\boldsymbol{A x}=\boldsymbol{b}$ 可化为 $\boldsymbol{P}^{\mathrm{T}} \boldsymbol{A P}\left(\boldsymbol{P}^{\mathrm{T}} \boldsymbol{x}\right)=\boldsymbol{P}^{\mathrm{T}} \boldsymbol{b}$，且记

$$\boldsymbol{y}=\boldsymbol{P}^{\mathrm{T}} \boldsymbol{x}=\binom{\boldsymbol{y}_{1}}{\boldsymbol{y}_{2}},\quad \boldsymbol{P}^{\mathrm{T}} \boldsymbol{b}=\binom{\boldsymbol{d}_{1}}{\boldsymbol{d}_{2}}$$

其中 $\boldsymbol{y}_{1},\boldsymbol{d}_{1}$ 为 $r$ 维向量。于是，求解 $\boldsymbol{A x}=\boldsymbol{b}$ 化为求解

$$\left\{\begin{aligned} \boldsymbol{A}_{11} \boldsymbol{y}_{1}+\boldsymbol{A}_{12} \boldsymbol{y}_{2} & =\boldsymbol{d}_{1},\\ \boldsymbol{A}_{22} \boldsymbol{y}_{2} & =\boldsymbol{d}_{2} \end{aligned}\right.$$

逐次超松弛迭代法

逐次超松弛迭代法（successive over relaxation method，简称 SOR 方法）是 Gauss-Seidel 方法的一种加速方法，是解大型稀疏矩阵方程组的有效方法之一，它具有计算公式简单，程序设计容易，占用计算机内存较少等优点，但需要选择好的加速因子（即最佳松弛因子）

设有方程组

$$\boldsymbol{A} \boldsymbol{x}=\boldsymbol{b}$$

其中 $\boldsymbol{A} \in \mathbf{R}^{n \times n}$ 为非奇异矩阵，且设 $a_{i i} \neq 0(i=1,2,\cdots,n)$，分解 $\boldsymbol{A}$ 为

$$A=D-L-U$$

设已知第 $k$ 次迭代向量 $\boldsymbol{x}^{(k)}$，及第 $k+1$ 次迭代向量 $\boldsymbol{x}^{(k+1)}$ 的分量 $x_{j}^{(k+1)} \quad(j=1,2,\cdots,i-1)$，要求计算分量 $x_{i}^{(k+1)}$。

首先用 Gauss-Seidel 迭代法定义辅助量

$$\tilde{x}_{i}^{(k+1)}=\frac{1}{a_{i i}}\left(b_{i}-\sum_{j=1}^{i-1} a_{i j} x_{j}^{(k+1)}-\sum_{j=i+1}^{n} a_{i j} x_{j}^{(k)}\right) \quad(i=1,2,\cdots,n)$$

再把 $x_{i}^{(k+1)}$ 取为 $x_{i}^{(k)}$ 与 $\tilde{x}_{i}^{(k+1)}$ 某个平均值（即加权平均），即

$$x_{i}^{(k+1)}=(1-\omega) x_{i}^{(k)}+\omega \tilde{x}_{i}^{(k+1)}=x_{i}^{(k)}+\omega\left(\tilde{x}_{i}^{(k+1)}-x_{i}^{(k)}\right)$$

代入即得到解方程组 $\boldsymbol{A x}=\boldsymbol{b}$ 的逐次超松弛迭代公式

$$\left\{\begin{array}{l} x_{i}^{(k+1)}=x_{i}^{(k)}+\frac{\omega}{a_{i i}}\left(b_{i}-\sum_{j=1}^{i-1} a_{i j} x_{j}^{(k+1)}-\sum_{j=i}^{n} a_{i j} x_{j}^{(k)}\right),\\ x^{(k)}=\left(x_{1}^{(k)},x_{2}^{(k)},\cdots,x_{n}^{(k)}\right)^{\mathrm{T}} \quad(k=0,1,\cdots;i=1,2,\cdots,n), \end{array}\right.$$

或写为

$$\left\{\begin{array}{l} x_{i}^{(k+1)}=x_{i}^{(k)}+\Delta x_{i} \quad(k=0,1,\cdots;i=1,2,\cdots,n)\\ \Delta x_{i}=\frac{\omega}{a_{i i}}\left(b_{i}-\sum_{j=1}^{i-1} a_{i j} x_{j}^{(k+1)}-\sum_{j=i}^{n} a_{i j} x_{j}^{(k)}\right) \end{array}\right.$$

矩阵形式

$$\boldsymbol{x}^{(k+1)}=\boldsymbol{L}_{\omega} \boldsymbol{x}^{(k)}+\boldsymbol{f}$$

其中

$$\boldsymbol{L}_{\omega}=(\boldsymbol{D}-\omega \boldsymbol{L})^{-1}[(1-\omega) \boldsymbol{D}+\omega \boldsymbol{U}],\quad \boldsymbol{f}=\omega(\boldsymbol{D}-\omega \boldsymbol{L})^{-1} \boldsymbol{b}$$

• 其中 $\omega$ 称为松弛因子
  • 当 $\omega=1$ 时，SOR 方法就是 Gauss-Seidel 迭代法
  • 当 $\omega<1$ 时，称为低松弛法
  • 当 $\omega>1$ 时，称为超松弛法
• 矩阵 $\boldsymbol{L}_{\omega}$ 称为 $SOR$ 方法的迭代矩阵

对原方程组应用 SOR 方法相当于对方程组 $\boldsymbol{x}=\boldsymbol{L}_{\omega} \boldsymbol{x}+\boldsymbol{f}$ 应用一般迭代法。

收敛性分析

设有线性方程组 $\boldsymbol{A} \boldsymbol{x}=\boldsymbol{b}$，且 $a_{i i} \neq 0(i=1,2,\cdots,n)$，则解方程组的 SOR 方法收敛的充要条件是

$$\rho\left(\boldsymbol{L}_{\omega}\right)<1$$

• SOR 方法收敛，则 $0<\omega<2$
• 如果 $\boldsymbol{A}$ 为对称正定矩阵，且 $0<\omega<2$，则 SOR 方法收敛

最佳松弛因子公式

$$\omega_{\mathrm{opt}}=\frac{2}{1+\sqrt{1-\rho^{2}\left(\boldsymbol{B}_{0}\right)}}$$