概论 误差

误差的来源

模型误差

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

观测误差

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

截断误差

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

• 例如，当函数 $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$，则可以确定一个缩小了的有根区间 $\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}$

计算步骤：

1. 准备待分解的矩阵 $\boldsymbol{A}$，始终在同一个矩阵里进行操作
2. 每次选择当前的行和列（直角形），划分成上面一行和左下一列
   
3. 首先计算行的数据（蓝色），等于原始数据减去当前行左侧的数据依次乘以当前列上面的数据的乘累加（注意如果数据不够则以少的算，第一次没数据，就不用做乘累加）
4. 再计算列的数据（橙色），等于原始数据减去当前行左侧的数据依次乘以当前列上面的数据的乘累加，最后再除以直角形的直角所在位置的元素
   
   
   
   
5. 最后，分拆出 $\boldsymbol{L}$ 和 $\boldsymbol{U}$ $$L=\left[\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 3 / 2 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 \\ 2 & 2 & 19 / 5 & 1 \end{array}\right] \quad U=\left[\begin{array}{cccc} 2 & 4 & 4 & 2 \\ 0 & -3 & 6 & 3 \\ 0 & 0 & -5 & 0 \\ 0 & 0 & 0 & -9 \end{array}\right]$$

Gauss 改写的三角分解法（适合手算）

1. 令 $\boldsymbol{A}$ 为待分解的矩阵，$\boldsymbol{L}$ 和 $\boldsymbol{P}$ 的初始值为单位阵，$\boldsymbol{U}$ 的初始值为 $\boldsymbol{A}$
2. 从第一行开始对 $\boldsymbol{U}$ 执行高斯消元
3. 若当前第 $k$ 行 $\boldsymbol{U}$ 的主元（第 $k$ 列）是 $0$，则需要交换第 $k$ 行和 $k+1$ 行，交换后 $\boldsymbol{L}$ 对角线左下角部分也需要交换第 $k$ 行和 $k+1$ 行（事实上此时只有第 $1 \sim k-1$ 列有非零值），然后交换 $\boldsymbol{P}$ 的第 $k$ 行和 $k+1$ 行
4. 确保主元非零后，将第 $k$ 行的若干倍减到下面的行使得当前主元前面全部变成 $0$，设把第 $k$ 行的 $\lambda$ 倍减到第 $m$ 行，则 $\boldsymbol{L}$ 的第 $m$ 行第 $k$ 列变成 $\lambda$
5. 完成高斯消元后，得到 $\boldsymbol{PA}=\boldsymbol{LU}$

不选主元的三角分解法（Doolittle 分解）

设 $\boldsymbol{A}$ 为非奇异矩阵，且有分解式 $\boldsymbol{A}=\boldsymbol{L} \boldsymbol{U}$（这要求 $\boldsymbol{A}$ 的对角线元素全部非零），其中 $\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_{r k} u_{k i}+u_{r i} \quad\left(\mathrm{当} r<k,l_{r k}=0 \mathrm{时}\right)$$

故

$$u_{r i}=a_{r i}-\sum_{k=1}^{r-1} l_{rk} 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}=\boldsymbol{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_{r i}=a_{r 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{P}$ 为单位矩阵应用和 $\boldsymbol{A}$ 一样的交换得到的矩阵。最后分解完成的式子为 $\boldsymbol{P A}=\boldsymbol{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}$ 的 Frobenius 范数：$\boldsymbol{F}(\boldsymbol{A})=\|\boldsymbol{A}\|_{F}=\left(\sum_{i, j=1}^{n} \boldsymbol{a}_{i j}^{2}\right)^{\frac{1}{2}}$

谱半径

设 $\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}$ 为

$$\boldsymbol{A}=\boldsymbol{D}-\boldsymbol{L}-\boldsymbol{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)}}$$

第二章 插值法

基础定义

设函数 $y=f(x)$ 在区间 $[a,b]$ 上有定义，且已知在点 $a \leqslant x_{0} \leqslant x_{1}<\cdots <x_{n} \leqslant b$ 上的值 $y_{0},y_{1},\cdots,y_{n}$，若存在一简单函数 $P(x)$，使

$$P\left(x_{i}\right)=y_{i} \quad(i=0,1,\cdots,n)$$

成立，就称 $P(x)$ 为 $f(x)$ 的插值函数，点 $x_{0},x_{1},\cdots,x_{n}$ 称为插值节点，包含插值节点的区间 $[a,b]$ 称为插值区间，求插值函数 $P(x)$ 的方法称为插值法。若 $P(x)$ 是次数不超过 $n$ 的代数多项式，即

$$P(x)=a_{0}+a_{1} x+\cdots+a_{n} x^{n},$$

其中 $a_{i}$ 为实数，就称 $P(x)$ 为插值多项式，相应的插值法称为多项式插值；若 $P(x)$ 为分段的多项式，就称之为分段插值；若 $P(x)$ 为三角多项式，就称之为三角插值

插值多项式的存在唯一性

设 $P(x)$ 是插值多项式，用 $H_{n}$ 代表所有次数不超过 $n$ 的多项式集合，于是 $P(x) \in H_{n}$。所谓插值多项式 $P(x)$ 存在唯一，就是指在集合 $H_{n}$ 中有且只有一个 $P(x)$ 满足式

$$P\left(x_{i}\right)=y_{i} \quad(i=0,1,\cdots,n)$$

展开得

$$\left\{\begin{array}{c} a_{0}+a_{1} x_{0}+\cdots+a_{n} x_{0}^{n}=y_{0},\\ a_{0}+a_{1} x_{1}+\cdots+a_{n} x_{1}^{n}=y_{1},\\ \vdots \\ a_{0}+a_{1} x_{n}+\cdots+a_{n} x_{n}^{n}=y_{n} \end{array}\right.$$

这是一个关于 $a_{0},a_{1},\cdots,a_{n}$ 的 $n+1$ 元线性方程组

要证明插值多项式的存在唯一，只要证明方程组的解存在唯一，也就是证明系数行列式

$$V_{n}\left(x_{0},x_{1},\cdots,x_{n}\right)=\left|\begin{array}{ccccc} 1 & x_{0} & x_{0}^{2} & \cdots & x_{0}^{n} \\ 1 & x_{1} & x_{1}^{2} & \cdots & x_{1}^{n} \\ \vdots & \vdots & \vdots & & \vdots \\ 1 & x_{n} & x_{n}^{2} & \cdots & x_{n}^{n} \end{array}\right|$$

不为零，其中 $V_{n}\left(x_{0},x_{1},\cdots,x_{n}\right)$ 称为 Vandermonde 行列式。利用行列式性质可得

$$V_{n}\left(x_{0},x_{1},\cdots,x_{n}\right)=\prod_{i=1}^{n} \prod_{j=0}^{i-1}\left(x_{i}-x_{j}\right)$$

由于 $i \neq j$ 时 $x_{i} \neq x_{j}$，上式所有因子 $x_{i}-x_{j} \neq 0$，于是 $V_{n}\left(x_{0},x_{1},\cdots,x_{n}\right) \neq 0$，故方程组存在唯一的一组解 $a_{0},a_{1},\cdots,a_{n}$。

线性插值与抛物插值

$$L_{1}(x)=y_{k}+\frac{y_{k+1}-y_{k}}{x_{k+1}-x_{k}}\left(x-x_{k}\right) \mathrm{(点斜式)}$$ $$L_{1}(x)=\frac{x_{k+1}-x}{x_{k+1}-x_{k}} y_{k}+\frac{x-x_{k}}{x_{k+1}-x_{k}} y_{k+1}\mathrm{(两点式)}$$ $$\begin{aligned} L_{2}(x) &=y_{k-1} \frac{\left(x-x_{k}\right)\left(x-x_{k+1}\right)}{\left(x_{k-1}-x_{k}\right)\left(x_{k-1}-x_{k+1}\right)}+y_{k} \frac{\left(x-x_{k-1}\right)\left(x-x_{k+1}\right)}{\left(x_{k}-x_{k-1}\right)\left(x_{k}-x_{k+1}\right)} \\ &+y_{k+1} \frac{\left(x-x_{k-1}\right)\left(x-x_{k}\right)}{\left(x_{k+1}-x_{k-1}\right)\left(x_{k+1}-x_{k}\right)} \end{aligned}$$

插值基函数

若 $n$ 次多项式 $l_{j}(x)(j=0,1,\cdots,n)$ 在 $n+1$ 个节点 $x_{0}<x_{1}<\cdots<x_{n}$ 上满足条件

$$l_{j}\left(x_{k}\right)=\left\{\begin{array}{ll} 1,& k=j,\\ 0,& k \neq j \end{array} \quad(j,k=0,1,\cdots,n),\right.$$

就称这 $n+1$ 个 $n$ 次多项式 $l_{0}(x),l_{1}(x),\cdots,l_{n}(x)$ 为节点 $x_{0},x_{1},\cdots,x_{n}$ 上的 $n$ 次插值基函数

• $n$ 次插值基函数为 $$l_{k}(x)=\frac{\left(x-x_{0}\right) \cdots\left(x-x_{k-1}\right)\left(x-x_{k+1}\right) \cdots\left(x-x_{n}\right)}{\left(x_{k}-x_{0}\right) \cdots\left(x_{k}-x_{k-1}\right)\left(x_{k}-x_{k+1}\right) \cdots\left(x_{k}-x_{n}\right)} \quad(k=0,1,\cdots,n)$$

拉格朗日插值

拉格朗日（Lagrange）插值多项式 $L_{n}(x)$ 可表示为

$$L_{n}(x)=\sum_{k=0}^{n} y_{k} l_{k}(x)$$

由 $l_{k}(x)$ 的定义知

$$L_{n}\left(x_{j}\right)=\sum_{k=0}^{n} y_{k} l_{k}\left(x_{j}\right)=y_{j} \quad(j=0,1,\cdots,n)$$

若引入记号

$$\omega_{n+1}(x)=\left(x-x_{0}\right)\left(x-x_{1}\right) \cdots\left(x-x_{n}\right),$$

容易求得

$$\omega_{n+1}^{\prime}\left(x_{k}\right)=\left(x_{k}-x_{0}\right) \cdots\left(x_{k}-x_{k-1}\right)\left(x_{k}-x_{k+1}\right) \cdots\left(x_{k}-x_{n}\right),$$

于是插值多项式也可表示成

$$L_{n}(x)=\sum_{k=0}^{n} y_{k} \frac{\omega_{n+1}(x)}{\left(x-x_{k}\right) \omega_{n+1}^{\prime}\left(x_{k}\right)}$$

• $n$ 次插值多项式 $L_{n}(x)$ 通常是次数为 $n$ 的多项式，特殊情况次数可能小于 $n$（例如，三点共线）

截断误差

若在 $[a,b]$ 上用 $L_{n}(x)$ 近似 $f(x)$，则其截断误差为 $R_{n}(x)=f(x)-L_{n}(x)$，$R_{n}(x)$ 也称为插值多项式的余项或插值余项

设

• $f^{(n)}(x)$ 在 $[a,b]$ 上连续；
• $f^{(n+1)}(x)$ 在 $(a,b)$ 内存在；
• 节点互不相同且 $a \leqslant x_{i} \leqslant b,\ i=0,1,\cdots,n$，

则 $\forall x \in[a,b]$，$\exists \xi \in(a,b)$：

$$R_{n}(x)=f(x)-L_{n}(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!} \omega_{n+1}(x)$$

• 当 $n=1$ 时，线性插值余项为 $$R_{1}(x)=\frac{1}{2} f^{\prime \prime}(\xi) \omega_{2}(x)=\frac{1}{2} f^{\prime \prime}(\xi)\left(x-x_{0}\right)\left(x-x_{1}\right),\quad \xi \in\left[x_{0},x_{1}\right]$$
• 当 $n=2$ 时，抛物插值的余项为 $$R_{2}(x)=\frac{1}{6} f^{\prime \prime \prime}(\xi)\left(x-x_{0}\right)\left(x-x_{1}\right)\left(x-x_{2}\right),\quad \xi \in\left[x_{0},x_{2}\right]$$
• 插值多项式 $L_{n}(x)$ 逼近 $f(x)$ 的截断误差限是 $$\left|R_{n}(x)\right| \leqslant \frac{M_{n+1}}{(n+1)!}\left|\omega_{n+1}(x)\right|$$ 其中 $$M_{n+1} = \max _{a \leqslant x \leqslant b}\left|f^{(n+1)}(x)\right|$$

差商

称 $f\left[x_{0},x_{k}\right]=\frac{f\left(x_{k}\right)-f\left(x_{0}\right)}{x_{k}-x_{0}}$ 为函数 $f(x)$ 关于点 $x_{0},x_{k}$ 的一阶差商，称

$$f\left[x_{0},x_{1},x_{k}\right]=\frac{f\left[x_{0},x_{k}\right]-f\left[x_{0},x_{1}\right]}{x_{k}-x_{1}}$$

为 $f(x)$ 关于点 $x_{0},x_{1},x_{k}$ 的二阶差商

一般地，称

$$f\left[x_{0},x_{1},\cdots,x_{k}\right]=\frac{f\left[x_{0},x_{1},\cdots,x_{k-2},x_{k}\right]-f\left[x_{0},x_{1},\cdots,x_{k-1}\right]}{x_{k}-x_{k-1}}$$

为 $f(x)$ 的 k 阶差商

性质

• k 阶差商可表示为函数值 $f\left(x_{0}\right),f\left(x_{1}\right),\cdots,f\left(x_{k}\right)$ 的线性组合，即 $$f\left[x_{0},x_{1},\cdots,x_{k}\right]=\sum_{j=0}^{k} \frac{f\left(x_{j}\right)}{\left(x_{j}-x_{0}\right) \cdots\left(x_{j}-x_{j-1}\right)\left(x_{j}-x_{j+1}\right) \cdots\left(x_{j}-x_{k}\right)}$$
• 差商与节点的排列次序无关，称为差商的对称性，即 $$f\left[x_{0},x_{1},\cdots,x_{k}\right]=f\left[x_{1},x_{0},x_{2},\cdots,x_{k}\right]=\cdots=f\left[x_{1},x_{2},\cdots,x_{k},x_{0}\right]$$
• $$f\left[x_{0},x_{1},\cdots,x_{k}\right]=\frac{f\left[x_{1},x_{2},\cdots,x_{k}\right]-f\left[x_{0},x_{1},\cdots,x_{k-1}\right]}{x_{k}-x_{0}}$$
• 若 $f(x)$ 在 $[a,b]$ 上存在 $n$ 阶导数，且节点 $x_{0},x_{1},\cdots,x_{n} \in[a,b]$，则 $n$ 阶差商与导数关系为 $$f\left[x_{0},x_{1},\cdots,x_{n}\right]=\frac{f^{(n)}(\xi)}{n!},\quad \xi \in[a,b]$$
• 若 $F(x)=c f(x)$，则 $$F\left[x_{0},x_{1},\cdots,x_{n}\right]=c f\left[x_{0},x_{1},\cdots,x_{n}\right]$$
• 若 $F(x)=f(x)+g(x)$，则 $$F\left[x_{0},x_{1},\cdots,x_{n}\right]=f\left[x_{0},x_{1},\cdots,x_{n}\right]+g\left[x_{0},x_{1},\cdots,x_{n}\right]$$

Newton 插值公式

$$f(x) = N_{n}(x) + R_{n}(x)$$

其中插值项

$$\begin{aligned} N_{n}(x)= & f\left(x_{0}\right)+f\left[x_{0},x_{1}\right]\left(x-x_{0}\right)+f\left[x_{0},x_{1},x_{2}\right]\left(x-x_{0}\right)\left(x-x_{1}\right)+\cdots \\ & +f\left[x_{0},x_{1},\cdots,x_{n}\right]\left(x-x_{0}\right)\left(x-x_{1}\right) \cdots\left(x-x_{n-1}\right) \end{aligned}$$

余项

$$R_{n}(x)= f(x)-N_{n}(x)=f\left[x,x_{0},\cdots,x_{n}\right] \omega_{n+1}(x)$$

分段低次插值

背景

• 拉格朗日插值多项式 $L_{n}(x)$ 近似 $f(x)$ 时，$L_{n}(x)$ 的次数不是越高精度越好，事实上当 $n \rightarrow \infty$ 时，$L_{n}(x)$ 不一定收敛到 $f(x)$
• 分段线性插值，就是将插值点用折线段连接起来逼近 $f(x)$

定义

设已知节点 $a=x_{0}<x_{1}<\cdots<x_{n}=b$ 上的函数值 $f_{0},f_{1},\cdots,f_{n}$，记 $h_{k}=x_{k+1}-x_{k}$，$h=\max _{k} \{h_{k}\}$ 称 $I_{h}(x)$ 为分段线性插值函数，如果满足：

• 记 $I_{h}(x) \in C[a,b]$；
• $I_{h}\left(x_{k}\right)=f_{k}(k=0,1,\cdots,n)$；
• $I_{h}(x)$ 在每个区间 $\left[x_{k},x_{k+1}\right]$ 上是线性函数，

则由定义可知，$I_{h}(x)$ 在每个小区间 $\left[x_{k},x_{k+1}\right]$ 上可表示为

$$I_{h}(x)=\frac{x-x_{k+1}}{x_{k}-x_{k+1}} f_{k}+\frac{x-x_{k}}{x_{k+1}-x_{k}} f_{k+1} \quad\left(x_{k} \leqslant x \leqslant x_{k+1}\right)$$

若用插值基函数表示，则 $I_{h}(x)$ 在整个区间 $[a,b]$ 上可表示为

$$I_{h}(x)=\sum_{j=0}^{n} f_{j} l_{j}(x)$$

其中基函数 $l_{j}(x)$ 满足条件 $l_{j}\left(x_{k}\right)=\delta_{j k}(j,k=0,1,\cdots,n)$，其形式是

$$l_{j}(x)=\left\{\begin{array}{ll} \frac{x-x_{j-1}}{x_{j}-x_{j-1}},& x_{j-1} \leqslant x \leqslant x_{j}\quad(j=0 \mathrm{略去}),\\ \frac{x_{j}-x_{j+1}}{x-x_{j+1}},& x_{j} \leqslant x \leqslant x_{j+1} \quad(j=n \mathrm{略去}),\\ 0,& x \in[a,b],x \notin\left[x_{j-1},x_{j+1}\right] \end{array}\right.$$

分段线性插值基函数 $l_{j}(x)$ 只在 $x_{j}$ 附近不为零，在其他地方均为零，这种性质称为局部非零性质

Hermite 插值

设在节点 $a \leqslant x_{0}<x_{1}<\cdots<x_{n} \leqslant b$ 上，$y_{j} =f\left(x_{j}\right)$，$m_{j}=f^{\prime}\left(x_{j}\right)(j=0,1,\cdots,n)$，要求插值多项式 $H(x)$，满足条件

$$H\left(x_{j}\right)=y_{j},\quad H^{\prime}\left(x_{j}\right)=m_{j} \quad(j=0,1,\cdots,n)$$

这里给出的 $2 n+2$ 个条件，可唯一确定一个次数不超过 $2 n+1$ 的多项式

$$H_{2 n+1}(x)=H(x)$$

其形式为

$$H_{2 n+1}(x)=a_{0}+a_{1} x+\cdots+a_{2 n+1} x^{2 n+1}$$

用求 Lagrange 插值多项式的基函数方法：先求插值基函数 $\alpha_{j}(x)$ 及 $\beta_{j}(x)(j=0,1,\cdots,n)$，共有 $2 n+2$ 个，每一个基函数都是 $2 n+1$ 次多项式，且满足条件

$$\left\{\begin{array}{l} \alpha_{j}\left(x_{k}\right)=\delta_{j k}=\left\{\begin{array}{ll} 0,& j \neq k,\\ 1,& j=k, \end{array} \quad \alpha_{j}^{\prime}\left(x_{k}\right)=0 \quad(j,k=0,1,\cdots,n),\right. \\ \beta_{j}\left(x_{k}\right)=0,\quad \beta_{j}^{\prime}\left(x_{k}\right)=\delta_{j k} \end{array}\right.$$

解得

$$\alpha_{j}(x)=\left[1-2\left(x-x_{j}\right) \sum_{\substack{k=0 \\ k \neq j}}^{n} \frac{1}{x_{j}-x_{k}}\right] l_{j}^{2}(x)$$ $$\beta_{j}(x)=\left(x-x_{j}\right) l_{j}^{2}(x)$$

于是插值多项式 $H(x)=H_{2 n+1}(x)$ 可写成用插值基函数表示的形式：

$$H_{2 n+1}(x)=\sum_{j=0}^{n}\left[y_{j} \alpha_{j}(x)+m_{j} \beta_{j}(x)\right]$$

重要特例是 $n=1$ 的情形，这时可取节点为 $x_{k}$ 及 $x_{k+1}$，插值多项式为 $H_{3}(x)$，满足条件

$$\left\{\begin{array}{l} H_{3}\left(x_{k}\right)=y_{k},\quad H_{3}\left(x_{k+1}\right)=y_{k+1},\\ H_{3}^{\prime}\left(x_{k}\right)=m_{k},\quad H_{3}^{\prime}\left(x_{k+1}\right)=m_{k+1} \end{array}\right.$$

基函数

$$\left\{\begin{array}{l} \alpha_{k}(x)=\left(1+2 \frac{x-x_{k}}{x_{k+1}-x_{k}}\right)\left(\frac{x-x_{k+1}}{x_{k}-x_{k+1}}\right)^{2} \\ \alpha_{k+1}(x)=\left(1+2 \frac{x-x_{k+1}}{x_{k}-x_{k+1}}\right)\left(\frac{x-x_{k}}{x_{k+1}-x_{k}}\right)^{2} \end{array}\right.$$ $$\left\{\begin{array}{l} \beta_{k}(x)=\left(x-x_{k}\right)\left(\frac{x-x_{k+1}}{x_{k}-x_{k+1}}\right)^{2} \\ \beta_{k+1}(x)=\left(x-x_{k+1}\right)\left(\frac{x-x_{k}}{x_{k+1}-x_{k}}\right)^{2} \end{array}\right.$$

插值多项式是

$$H_{3}(x)=y_{k} \alpha_{k}(x)+y_{k+1} \alpha_{k+1}(x)+m_{k} \beta_{k}(x)+m_{k+1} \beta_{k+1}(x)$$

其余项 $R_{3}(x)=f(x)-H_{3}(x)$

$$R_{3}(x)=\frac{1}{4!} f^{(4)}(\xi)\left(x-x_{k}\right)^{2}\left(x-x_{k+1}\right)^{2}$$

分段三次 Hermite 插值

分段线性插值函数 $I_{h}(x)$ 的导数是间断的，若在节点 $x_{k}(k=0,1,\cdots,n)$ 上除已知函数值 $f_{k}$ 外还给出导数值 $f_{k}^{\prime}=m_{k}(k=0,1,\cdots,n)$，就可构造一个导数连续的分段插值函数 $I_{h}(x)$，它满足如下条件：

• $I_{h}(x) \in C^{1}[a,b]$（$C^{1}[a,b]$ 代表区间 $[a,b]$ 上一阶导数连续的函数集合）
• $I_{h}\left(x_{k}\right)=f_{k},I_{h}^{\prime}\left(x_{k}\right)=f_{k}^{\prime} \quad(k=0,1,\cdots,n)$
• $I_{h}(x)$ 在每个小区间 $\left[x_{k},x_{k+1}\right]$ 上是三次多项式 $$\alpha_{j}(x)= \left\{\begin{array}{ll} \left(\frac{x-x_{j-1}}{x_{j}-x_{j-1}}\right)^{2}\left(1+2 \frac{x-x_{j}}{x_{j-1}-x_{j}}\right),& x_{j-1} \leqslant x \leqslant x_{j} \quad(j=0 \mathrm{略去}),\\ \left(\frac{x-x_{j+1}}{x_{j}-x_{j+1}}\right)^{2}\left(1+2 \frac{x-x_{j}}{x_{j+1}-x_{j}}\right),& x_{j} \leqslant x \leqslant x_{j+1} \quad(j=n \mathrm{略去}),\\ 0,& \mathrm{其他} \end{array}\right.$$ $$\beta_{j}(x)= \left\{\begin{array}{ll} \left(\frac{x-x_{j-1}}{x_{j}-x_{j-1}}\right)^{2}\left(x-x_{j}\right),& x_{j-1} \leqslant x \leqslant x_{j} \quad(j=0 \mathrm{略去}),\\ \left(\frac{x-x_{j+1}}{x_{j}-x_{j+1}}\right)^{2}\left(x-x_{j}\right),& x_{j} \leqslant x \leqslant x_{j+1} \quad(j=n \mathrm{略去}),\\ 0,& \mathrm{其他} \end{array}\right.$$ 于是 $I_{h}(x)$ 可表示为 $$I_{h}(x)=\sum_{j=0}^{n}\left[f_{j} \alpha_{j}(x)+f_{j}^{\prime} \beta_{j}(x)\right]$$ 由于 $\alpha_{j}(x),\beta_{j}(x)$ 的局部非零性质，当 $x \in\left[x_{k},x_{k+1}\right]$ 时，只有 $\alpha_{k}(x),\alpha_{k+1}(x)$，$\beta_{k}(x),\beta_{k+1}(x)$ 不为零，于是 $I_{h}(x)$ 也可表示为 $$I_{h}(x)=f_{k} \alpha_{k}(x)+f_{k+1} \alpha_{k+1}(x)+f_{k}^{\prime} \beta_{k}(x)+f_{k+1}^{\prime} \beta_{k+1}(x) \quad\left(x_{k} \leqslant x \leqslant x_{k+1}\right)$$

三次样条曲线

定义

若函数 $S(x) \in C^{2}[a,b]$，且在每个小区间 $\left[x_{j},x_{j+1}\right]$ 上是三次多项式，其中 $a=x_{0}<x_{1}<\cdots<x_{n}=b$ 是给定节点，则称 $S(x)$ 是节点 $x_{0},x_{1},\cdots,x_{n}$ 上的三次样条函数。若在节点 $x_{j}$ 上给定函数值 $y_{j}=f\left(x_{j}\right)(j=0,1,\cdots,n)$，且

$$S\left(x_{j}\right)=y_{j} \quad(j=0,1,\cdots,n)$$

成立，则称 $S(x)$ 为三次样条插值函数。

解析条件

从定义知，要求出 $S(x)$，在每个小区间 $\left[x_{j},x_{j+1}\right]$ 上要确定 $4$ 个待定系数，共有 $n$ 个小区间，故应确定 $4 n$ 个参数。根据 $S(x)$ 在 $[a,b]$ 上二阶导数连续，在节点 $x_{j} (j=1,2,\cdots,n-1)$ 处应满足连续性条件

$$\begin{aligned} S\left(x_{j}-0\right)= & S\left(x_{j}+0\right),\quad S^{\prime}\left(x_{j}-0\right)=S^{\prime}\left(x_{j}+0\right),\\ & S^{\prime \prime}\left(x_{j}-0\right)=S^{\prime \prime}\left(x_{j}+0\right), \end{aligned}$$

共有 $3 n-3$ 个条件，再加上 $S(x)$ 满足插值条件，共有 $4 n-2$ 个条件，因此还需要 $2$ 个条件才能确定 $S(x)$ 通常可在区间 $[a,b]$ 端点 $a=x_{0},b=x_{n}$ 上各加一个条件（称为边界条件）。

• 已知两端的一阶导数值，即 $$\left\{\begin{array}{l} S^{\prime}\left(x_{0}\right)=f_{0}^{\prime},\\ S^{\prime}\left(x_{n}\right)=f_{n}^{\prime} \end{array}\right.$$
• 两端的二阶导数已知，即 $$\left\{\begin{array}{l} S^{\prime \prime}\left(x_{0}\right)=f_{0}^{\prime \prime},\\ S^{\prime \prime}\left(x_{n}\right)=f_{n}^{\prime \prime} \end{array}\right.$$
• 自然边界条件 $$S^{\prime \prime}\left(x_{0}\right)=S^{\prime \prime}\left(x_{n}\right)=0$$
• 当 $f(x)$ 是以 $x_{n}-x_{0}$ 为周期的周期函数时，则要求 $S(x)$ 也是周期函数，这时边界条件应满足 $$\left\{\begin{array}{l} S\left(x_{0}+0\right)=S\left(x_{n}-0\right),\\ S^{\prime}\left(x_{0}+0\right)=S^{\prime}\left(x_{n}-0\right),\\ S^{\prime \prime}\left(x_{0}+0\right)=S^{\prime \prime}\left(x_{n}-0\right), \end{array}\right.$$ 这样确定的样条函数 $S(x)$ 称为周期样条函数

三转角方程

若假定 $S^{\prime}(x)$ 在节点 $x_{j}$ 处的值为 $S^{\prime}\left(x_{j}\right)=m_{j} \quad(j=0,1,\cdots,n)$，$h_{j}=x_{j+1}-x_{j}$，则由分段三次 Hermite 插值式可得目标函数

$$S(x)=\sum_{j=0}^{n}\left[y_{j} \alpha_{j}(x)+m_{j} \beta_{j}(x)\right]$$

应用已知条件，方程组可简化为

$$\lambda_{j} m_{j-1}+2 m_{j}+\mu_{j} m_{j+1}=g_{j} \quad(j=1,2,\cdots,n-1)$$

其中

$$\begin{array}{c} \lambda_{j}=\frac{h_{j}}{h_{j-1}+h_{j}},\quad \mu_{j}=\frac{h_{j-1}}{h_{j-1}+h_{j}} \quad(j=1,2,\cdots,n-1),\\ g_{j}=3\left(\lambda_{j} f\left[x_{j-1},x_{j}\right]+\mu_{j} f\left[x_{j},x_{j+1}\right]\right) \quad(j=1,2,\cdots,n-1) \end{array}$$

这个方程是关于 $n$ 个未知数 $m_{0},m_{1},\cdots,m_{n}$ 的 $n-1$ 个方程。

• 若已知两端的一阶导数值，则 $$\begin{array}{c} m_{0}=f_{0}^{\prime}\\ m_{n}=f_{n}^{\prime} \end{array}$$ 那么方程为只含 $m_{1},\cdots,m_{n-1}$ 的 $n-1$ 个方程，写成矩阵形式 $$\left(\begin{array}{cccccc} 2 & \mu_{1} & 0 & \cdots & \cdots & 0 \\ \lambda_{2} & 2 & \mu_{2} & \ddots & & \vdots \\ 0 & \lambda_{3} & 2 & \mu_{3} & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & 0 \\ 0 & & \ddots & \lambda_{n-2} & 2 & \mu_{n-2} \\ 0 & \cdots & \cdots & 0 & \lambda_{n-1} & 2 \end{array}\right)\left(\begin{array}{c} m_{1} \\ m_{2} \\ m_{3} \\ \vdots \\ m_{n-2} \\ m_{n-1} \end{array}\right)=\left(\begin{array}{c} g_{1}-\lambda_{1} f_{0}^{\prime} \\ g_{2} \\ g_{3} \\ \vdots \\ g_{n-2} \\ g_{n-1}-\mu_{n-1} f_{n}^{\prime} \end{array}\right)$$
• 若已知两端的二阶导数值，则 $$\left\{\begin{array}{l} 2 m_{0}+m_{1}=3 f\left[x_{0},x_{1}\right]-\frac{h_{0}}{2} f_{0}^{\prime \prime}=g_{0} \\ m_{n-1}+2 m_{n}=3 f\left[x_{n-1},x_{n}\right]+\frac{h_{n-1}}{2} f_{n}^{\prime \prime}=g_{n} \end{array}\right.$$ 于是，合并后用矩阵形式表示为 $$\left(\begin{array}{cccccc} 2 & 1 & 0 & \cdots & \cdots & 0 \\ \lambda_{1} & 2 & \mu_{1} & \ddots & & \vdots \\ 0 & \lambda_{2} & 2 & \mu_{2} & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & 0 \\ \vdots & & \ddots & \lambda_{n-1} & 2 & \mu_{n-1} \\ 0 & \cdots & \cdots & 0 & 1 & 2 \end{array}\right)\left(\begin{array}{c} m_{0} \\ m_{1} \\ m_{2} \\ \vdots \\ m_{n-1} \\ m_{n} \end{array}\right)=\left(\begin{array}{c} g_{0} \\ g_{1} \\ g_{2} \\ \vdots \\ g_{n-1} \\ g_{n} \end{array}\right)$$
• 若边界条件为周期样条，则 $$\begin{array}{c} m_{0}=m_{n} \\ \mu_{n} m_{1}+\lambda_{n} m_{n-1}+2 m_{n}=g_{n} \\ \end{array}$$ 其中 $$\begin{array}{c} \mu_{n}=\frac{h_{n-1}}{h_{0}+h_{n-1}},\quad \lambda_{n}=\frac{h_{n-1}}{h_{0}+h_{n-1}},\\ g_{n}=3\left(\mu_{n} f\left[x_{0},x_{1}\right]+\lambda_{n} f\left[x_{n-1},x_{n}\right]\right) \end{array}$$ 合并后的矩阵形式 $$\left(\begin{array}{cccccc} 2 & \mu_{1} & \cdots & 0 & \lambda_{1} \\ \lambda_{2} & 2 & \mu_{2} & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \ddots & \lambda_{n-1} & 2 & \mu_{n-1} \\ \mu_{n} & 0 & \cdots & \lambda_{n} & 2 \end{array}\right)\left(\begin{array}{c} m_{1} \\ m_{2} \\ \vdots \\ m_{n-1} \\ m_{n} \end{array}\right)=\left(\begin{array}{c} g_{1} \\ g_{2} \\ \vdots \\ g_{n-1} \\ g_{n} \end{array}\right)$$

收敛性

若 $f(x) \in C[a,b]$，$S(x)$ 是以 $a=x_{0}<x_{1}<\cdots<x_{n}=b$ 为节点，满足自然边界条件的三次样条插值函数，令

$$h_{j}=x_{j+1}-x_{j},\quad h=\max _{0 \leqslant j \leqslant n-1} h_{j},\quad \delta=\min _{0 \leqslant j \leqslant n-1} h_{j},$$

设 $h / \delta \leqslant c<\infty$，则 $S(x)$ 在 $[a,b]$ 上一致收敛到 $f(x)$

第三章 逼近

逼近的定义

对于函数类 $A$ 中给定的函数 $f(x)$，要求在另一类较简单的便于计算的函数类 $B$ 中，求函数 $P(x)\in B \subset A$，使 $P(x)$ 与 $f(x)$ 之差在某种度量意义下最小

• 函数类 $A$ 通常是区间 $[a,b]$ 上的连续函数,记作 $C[a,b]$
• 函数类 $B$ 通常是代数多项式、分式有理函数或三角多项式
• 度量标准最常用的有两种：
  • 一致逼近（均匀逼近） $$\|f(x)-P(x)\|_{\infty}=\max _{a \leqslant x \leqslant b}|f(x)-P(x)|$$
  • 均方逼近（平方逼近）
    在这种度量意义下的函数逼近称为；另一种度量标准是 $$\|f(x)-P(x)\|_{2}=\sqrt{\int_{a}^{b}[f(x)-P(x)]^{2} \mathrm{~d} x}$$

Weierstrass 定理

设 $f(x) \in C[a,b]$，则对于任何 $\varepsilon>0$，总存在一个代数多项式 $P(x)$，使

$$\|f(x)-P(x)\|_{\infty}<\varepsilon$$

在 $[a,b]$ 上一致成立

函数范数

对于定义在一个区间 $[a,b]$ 上的实连续函数 $f(x)\in C[a,b]$，$L^{p}$ 范数（$p \geq 1$）定义为：

$$\|f\|_{p}=\left(\int_{a}^{b}|f(x)|^{p} \mathrm{d} x\right)^{\frac{1}{p}}$$

• 当 $p \rightarrow \infty$ 时，定义 $\|f\|_{\infty}=\sup _{x \in[a,b]}|f(x)|$。

三种常用范数：

• $L^{1}$ 范数：$\|f\|_{1}=\int_{a}^{b}|f(x)| \mathrm{d}x$ 表示函数绝对值的积分。
• $L^{2}$ 范数：$\|f\|_{2}=\sqrt{\int_{a}^{b}|f(x)|^{2} \mathrm{d} x}$ 表示对应函数的欧几里得长度，经常用于内积空间。
• $L^{\infty}$ 范数：$\|f\|_{\infty}=\max _{x \in[a,b]}|f(x)|$ 表示函数在区间 $[a,b]$ 上的最大绝对值。

性质

• $\|f\| \geqslant 0$，当且仅当 $f \equiv 0$ 时才有 $\|f\|=0$
• $\|a f\|=|a|\|f\|$ 对于任意 $f \in C[a,b]$ 成立，$a$ 为任意实数
• 对于任意 $f,g \in C[a,b]$，有 $$\|f+g\| \leqslant\|f\|+\|g\|$$

函数的距离

与向量空间类似，当 $f,g \in C[a,b]$ 时，定义 $f$ 与 $g$ 的距离为

$$D(f,g)=\|f-g\|_{\infty}$$

于是可得到

$$\begin{array}{c} D(f,g) \leqslant D(f,h)+D(h,g) \\ \left|\|f\|_{\infty}-\|g\|_{\infty}\right| \leqslant\|f-g\|_{\infty} \end{array}$$

最佳一致逼近问题

逼近多项式

记次数不大于 $n$ 的多项式集合为 $H_{n}$，$H_{n} \subseteq C[a,b]$。又记 $H_{n}=\operatorname{span}\left\{1,x,\cdots,x^{n}\right\}$，其中 $1,x,\cdots,x^{n}$ 是 $[a,b]$ 上一组线性无关的函数组，是 $H_{n}$ 中的一组基，$H_{n}$ 中的元素 $P_{n}(x)$ 可表示为

$$P_{n}(x)=a_{0}+a_{1} x+\cdots+a_{n} x^{n}$$

偏差

设 $P_{n}(x) \in H_{n}$，$f(x) \in C[a,b]$，称

$$\Delta\left(f,P_{n}\right)=\left\|f-P_{n}\right\|_{\infty}=\max _{a \leqslant x \leqslant b}\left|f(x)-P_{n}(x)\right|$$

为 $f(x)$ 与 $P_{n}(x)$ 在 $[a,b]$ 上的偏差

偏差的全体组成一个集合，记作 $\left\{\Delta\left(f,P_{n}\right)\right\}$，下界是 $0$，下确界记为

$$E_{n}=\inf _{P_{n} \in H_{n}}\left\{\Delta\left(f,P_{n}\right)\right\}=\inf _{P_{n} \in H_{n}} \max _{a \leqslant x \leqslant b}\left|f(x)-P_{n}(x)\right|$$

则称 $E_{n}$ 为 $f(x)$ 在 $[a,b]$ 上的最小偏差，这个最小偏差在逼近多项式为最佳逼近多项式时取得

最佳一致逼近多项式

假定 $f(x) \in C[a,b]$，若存在 $P_{n}^{*}(x) \in H_{n}$ 使得

$$\Delta\left(f,P_{n}^{*}\right)=E_{n},$$

则称 $P_{n}^{*}(x)$ 是 $f(x)$ 在 $[a,b]$ 上的最佳一致逼近多项式或最小偏差逼近多项式，简称最佳逼近多项式

存在性：若 $f(x) \in C[a,b]$，则总存在 $P_{n}^{*}(x) \in H_{n}$，使

$$\left\|f(x)-P_{n}^{*}(x)\right\|_{\infty}=E_{n}$$

也就是最佳逼近多项式一定存在。

偏差点

设 $f(x) \in C[a,b]$，$P(x) \in H_{n}$，若在 $x=x_{0}$ 上有

$$\left|P\left(x_{0}\right)-f\left(x_{0}\right)\right|=\max _{a \leqslant x \leqslant b}|P(x)-f(x)|=\mu,$$

则称 $x_{0}$ 是 $P(x)$ 的偏差点

• 若 $P\left(x_{0}\right)-f\left(x_{0}\right)=\mu$，则称 $x_{0}$ 为正偏差点
• 若 $P\left(x_{0}\right)-f\left(x_{0}\right)=-\mu$，则称 $x_{0}$ 为负偏差点

存在性：由于函数 $P(x)-f(x)$ 在 $[a,b]$ 上连续，因此，至少存在一个点 $x_{0} \in[a,b]$，使得

$$\left|P\left(x_{0}\right)-f\left(x_{0}\right)\right|=\mu,$$

Chebyshev（切比雪夫）定理

$P(x) \in H_{n}$ 是 $f(x) \in C[a,b]$ 的最佳逼近多项式的充分必要条件是 $P(x)$ 在 $[a,b]$ 上至少有 $n+2$ 个轮流为正、负的偏差点，即有 $n+2$ 个点 $a \leqslant x_{1}< x_{2}<\cdots<x_{n+2} \leqslant b$，使

$$\begin{array}{l} P\left(x_{k}\right)-f\left(x_{k}\right)=(-1)^{k} \sigma\|P(x)-f(x)\|_{\infty},\\ \sigma= \pm 1,\quad k=1,2,\cdots,n+2 \end{array}$$

这样的点组称为 Chebyshev 交错点组

• 若 $f(x) \in C[a,b]$，则在 $H_{n}$ 中存在唯一的最佳逼近多项式
• 若 $f(x) \in C[a,b]$，则其最佳逼近多项式 $P_{n}^{*}(x) \in H_{n}$ 是 $f(x)$ 的一个 Lagrange 插值多项式。

内积空间

权函数

设在区间 $(a,b)$ 内，非负函数 $\rho(x)$ 满足以下条件，就称 $\rho(x)$ 为区间 $(a,b)$ 内的权函数：

• $\int_{a}^{b}|x|^{n} \rho(x) \mathrm{d} x(n=0,1,\cdots)$ 存在；
• 对于非负的连续函数 $g(x)$，若 $$\int_{a}^{b} g(x) \rho(x) \mathrm{d} x=0,$$ 则在 $(a,b)$ 内 $g(x) \equiv 0$

内积

设 $f(x),g(x) \in C[a,b]$，$\rho(x)$ 是 $[a,b]$ 上的权函数，积分

$$(f,g)=\int_{a}^{b} \rho(x) f(x) g(x) \mathrm{d} x$$

称为函数 $f(x)$ 与 $g(x)$ 在 $[a,b]$ 上的内积

内积满足下列四条公理：

• $(f,g)=(g,f)$；
• $(c f,g)=c(f,g)$，$c$ 为常数；
• $\left(f_{1}+f_{2},g\right)=\left(f_{1},g\right)+\left(f_{2},g\right)$；
• $(f,f) \geqslant 0$，当且仅当 $f=0$ 时 $(f,f)=0$

满足内积定义的函数空间称为内积空间

• 连续函数空间 $C[a,b]$ 上定义了内积就形成一个内积空间

欧几里得范数

$f(x) \in C[a,b]$，则

$$\|f\|_{2}=\sqrt{\int_{a}^{b} \rho(x) f^{2}(x) \mathrm{d} x}=\sqrt{(f,f)}$$

称为 $f(x)$ 的 Euclid 范数

对于任何 $f,g \in C[a,b]$，下列结论成立：

• $\|(f,g)\| \leqslant\|f\|_{2}\|g\|_{2}$ （Cauchy-Schwarz 不等式）
• $\|f+g\|_{2} \leqslant\|f\|_{2}+\|g\|_{2}$（三角不等式）
• $\|f+g\|_{2}^{2}+\|f-g\|_{2}^{2}=2\left(\|f\|_{2}^{2}+\|g\|_{2}^{2}\right)$（平行四边形定律）

正交函数

若 $f(x),g(x) \in C[a,b]$ 满足

$$(f,g)=\int_{a}^{b} \rho(x) f(x) g(x) \mathrm{d} x=0,$$

则称 $f$ 与 $g$ 在 $[a,b]$ 上带权 $\rho(x)$ 正交

若函数族 $\varphi_{0}(x),\varphi_{1}(x),\cdots,\varphi_{n}(x),\cdots$ 满足关系

$$\left(\varphi_{j},\varphi_{k}\right)=\int_{a}^{b} \rho(x) \varphi_{j}(x) \varphi_{k}(x) \mathrm{d} x=\left\{\begin{array}{ll} 0,& j \neq k,\\ A_{k}>0,& j=k, \end{array}\right.$$

则称 $\left\{\varphi_{k}\right\}$ 是 $[a,b]$ 上带权 $\rho(x)$ 的正交函数族

若 $A_{k} \equiv 1$，则称 $\left\{\varphi_{k}\right\}$ 为标准正交函数族

例如，三角函数族

$$1,\cos x,\sin x,\cos 2 x,\sin 2 x,\cdots$$

就是在区间 $[-\pi,\pi]$ 上的正交函数族（权 $\rho(x) \equiv 1$），其

$$\begin{array}{l} (1,1)=2 \pi,\\ (\sin n x,\sin m x)=\int_{-\pi}^{\pi} \sin n x \sin m x \mathrm{~d} x=\left\{\begin{array}{ll} \pi,& m=n,\\ 0,& m \neq n \end{array},(n,m=1,2,\cdots),\right. \\ (\cos n x,\cos m x)=\int_{-\pi}^{\pi} \cos n x \cos m x \mathrm{~d} x=\left\{\begin{array}{ll} \pi,& m=n,\\ 0,& m \neq n \end{array}(n,m=1,2,\cdots),\right. \\ (\cos n x,\sin m x)=\int_{-\pi}^{\pi} \cos n x \sin m x \mathrm{~d} x=0 \quad(n,m=1,2,\cdots) \end{array}$$

在 $\mathbf{R}^{n}$ 空间中，任一向量都可用它的一组线性无关的基表示，内积空间的任一元素 $f(x) \in C[a,b]$ 也同样可用线性无关的基表示

线性无关函数族

设 $\varphi_{0}(x),\varphi_{1}(x),\cdots,\varphi_{n-1}(x)$ 在 $[a,b]$ 上连续，如果

$$a_{0} \varphi_{0}(x)+a_{1} \varphi_{1}(x)+\cdots+a_{n-1} \varphi_{n-1}(x)=0$$

当且仅当 $a_{0}=a_{1}=\cdots=a_{n-1}=0$ 时成立，则称 $\varphi_{0},\varphi_{1},\cdots,\varphi_{n-1}$ 在 $[a,b]$ 上是线性无关的

若函数族 $\left\{\varphi_{k}\right\}(k=0,1,\cdots)$ 中的任何有限个 $\varphi_{k}$ 线性无关，则称 $\left\{\varphi_{k}\right\}$ 为线性无关函数族

若 $\varphi_{0}(x),\varphi_{1}(x),\cdots,\varphi_{n-1}(x)$ 是 $[a,b]$ 上的线性无关函数，且 $a_{0},a_{1},\cdots,a_{n-1}$ 是任意实数，则

$$S(x)=a_{0} \varphi_{0}(x)+a_{1} \varphi_{1}(x)+\cdots+a_{n-1} \varphi_{n-1}(x)$$

的全体是 $C[a,b]$ 中的一个子集，记作

$$\Phi=\operatorname{span}\left\{\varphi_{0},\varphi_{1},\cdots,\varphi_{n-1}\right\}$$

下面给出判断函数族 $\left\{\varphi_{k}\right\}(k=0,1,\cdots,n-1)$ 线性无关的充要条件

$\varphi_{0}(x),\varphi_{1}(x),\cdots,\varphi_{n-1}(x)$ 在 $[a,b]$ 上线性无关的充分必要条件是它的 Cramer 行列式 $G_{n-1} \neq 0$，其中

$$G_{n-1}=G\left(\varphi_{0},\varphi_{1},\cdots,\varphi_{n-1}\right)=\left|\begin{array}{cccc} \left(\varphi_{0},\varphi_{0}\right) & \left(\varphi_{0},\varphi_{1}\right) & \cdots & \left(\varphi_{0},\varphi_{n-1}\right) \\ \left(\varphi_{1},\varphi_{0}\right) & \left(\varphi_{1},\varphi_{1}\right) & \cdots & \left(\varphi_{1},\varphi_{n-1}\right) \\ \vdots & \vdots & & \vdots \\ \left(\varphi_{n-1},\varphi_{0}\right) & \left(\varphi_{n-1},\varphi_{1}\right) & \cdots & \left(\varphi_{n-1},\varphi_{n-1}\right) \end{array}\right|$$

正交多项式

设 $g_{n}(x)$ 是首项系数 $a_{n} \neq 0$ 的 $n$ 次多项式，如果多项式序列 $g_{0}(x),g_{1}(x),\cdots$ 满足

$$\left(g_{j},g_{k}\right)=\int_{a}^{b} \rho(x) g_{j}(x) g_{k}(x) \mathrm{d} x=\left\{\begin{array}{cc} 0,& j \neq k,\\ A_{k}>0,& j=k \end{array} \quad(j,k=0,1,\cdots),\right.$$

则称多项式序列 $g_{0}(x),g_{1}(x),\cdots$ 在 $[a,b]$ 上带权 $\rho(x)$ 正交，并称 $g_{n}(x)$ 是 $[a,b]$ 上带权 $\rho(x)$ 的 n 次正交多项式

一般来说，当区间 $[a,b]$ 及权函数 $\rho(x)$ 给定以后，可以由线性无关的一组基 $\left\{1,x,x^{2},\cdots,x^{n},\cdots\right\}$ 并利用正交化方法构造出正交多项式

$$g_{0}(x)=1,\quad g_{n}(x)=x^{n}-\sum_{k=0}^{n-1} \frac{\left(x^{n},g_{k}\right)}{\left(g_{k},g_{k}\right)} \cdot g_{k}(x),\quad k=1,2,\cdots$$

这样构造的正交多项式有以下性质：

• $g_{n}(x)$ 是最高项系数为 $1$ 的 $n$ 次多项式。
• 任一 $n$ 次多项式 $P_{n}(x) \in H_{n}$ 均可表示为 $g_{0}(x),g_{1}(x),\cdots,g_{n}(x)$ 的线性组合
• 当 $n \neq m$ 时，$\left(g_{n},g_{m}\right)=0$ 且 $g_{n}(x)$ 与任一次数小于 $n$ 的多项式正交
• 有递推关系 $$g_{n+1}(x)=\left(x-\alpha_{n}\right) g_{n}(x)-\beta_{n} g_{n-1}(x),\quad n=0,1,\cdots,$$ 其中 $$\begin{array}{c} g_{0}(x)=1, \quad g_{n-1}(x)=0, \\ \alpha_{n}=\frac{\left(x g_{n}, g_{n}\right)}{\left(g_{n}, g_{n}\right)}, \quad n=0,1, \cdots, \\ \beta_{n}=\frac{\left(g_{n},g_{n}\right)}{\left(g_{n-1},g_{n-1}\right)},\quad n=1,2,\cdots \end{array}$$ 这里 $\left(x g_{n},g_{n}\right)=\int_{a}^{b} x g_{n}^{2}(x) \rho(x) \mathrm{d} x$
• 设 $g_{0}(x),g_{1}(x),\cdots$ 是在 $[a,b]$ 上带权 $\rho(x)$ 的正交多项式序列，则 $g_{n}(x)(n \geqslant 1)$ 的 $n$ 个根都是单重实根，且都在区间 $(a,b)$ 内

Legendre（勒让德）多项式

当区间为 $[-1,1]$、权函数 $\rho(x) \equiv 1$ 时，由 $\left\{1,x,x^{2},\cdots,x^{n},\cdots\right\}$ 正交化得到的多项式就称为 Legendre 多项式，并用 $P_{0}(x),P_{1}(x),\cdots,P_{n}(x),\cdots$ 表示

$$P_{0}(x)=1,\quad P_{n}(x)=\frac{1}{2^{n} n!} \frac{\mathrm{d}^{n}}{\mathrm{~d} x^{n}}\left\{\left(x^{2}-1\right)^{n}\right\} \quad(n=1,2,\cdots)$$

由于 $\left(x^{2}-1\right)^{n}$ 是 $2 n$ 次多项式，求 $n$ 阶导数后得

$$P_{n}(x)=\frac{1}{2^{2} n!}(2 n)(2 n-1) \cdots(n+1) x^{n}+a_{n-1} x^{n-1}+\cdots+a_{0},$$

于是得首项 $x^{n}$ 的系数 $a_{n}=\frac{(2 n)!}{2^{n}(n!)^{2}}$，显然最高项系数为 $1$ 的 Legendre 多项式为

$$\widetilde{P}_{n}(x)=\frac{n!}{(2 n)!} \frac{\mathrm{d}^{n}}{\mathrm{~d} x^{n}}\left[\left(x^{2}-1\right)^{n}\right]$$

性质

• 正交性 $$\int_{-1}^{1} P_{n}(x) P_{m}(x) \mathrm{d} x=\left\{\begin{array}{cc} 0,& m \neq n,\\ \frac{2}{2 n+1},& m=n \end{array}\right.$$
• 奇偶性 $$P_{n}(-x)=(-1)^{n} P_{n}(x)$$
• 递推关系 $$(n+1) P_{n+1}(x)=(2 n+1) x P_{n}(x)-n P_{n-1}(x) \quad(n=1,2, \cdots)$$
• 在所有最高项系数为 $1$ 的 $n$ 次多项式中，Legendre 多项式 $\widetilde{P}_{n}(x)$ 在 $[-1,1]$ 上与零的平方误差最小
• $P_{n}(x)$ 在区间 $(-1,1)$ 内有 $n$ 个不同的实零点。

低次多项式

$$\begin{array}{l} P_{0}(x)=1 \\ P_{1}(x)=x \\ P_{2}(x)=\left(3 x^{2}-1\right) / 2 \\ P_{3}(x)=\left(5 x^{3}-3 x\right) / 2 \\ P_{4}(x)=\left(35 x^{4}-30 x^{2}+3\right) / 8 \\ P_{5}(x)=\left(63 x^{5}-70 x^{3}+15 x\right) / 8 \\ P_{6}(x)=\left(231 x^{6}-315 x^{4}+105 x^{2}-5\right) / 16 \end{array}$$

其中

• $P_{2}(x)$ 的根为 $\pm \frac{1}{\sqrt{3}}$
• $P_{3}(x)$ 的根为 $0$, $\pm \sqrt{\frac{3}{5}}$
• $P_{4}(x)$ 的根为 $\pm \sqrt{\frac{3}{7} \pm \frac{2}{7}\sqrt{\frac{6}{5}}}$
• $P_{5}(x)$ 的根为 $0$, $\pm \frac{1}{3}\sqrt{5 \pm 2\sqrt{\frac{10}{7}}}$

Chebyshev（切比雪夫）多项式

当权函数 $\rho(x)=\frac{1}{\sqrt{1-x^{2}}}$，区间为 $[-1,1]$ 时，由序列 $\left\{1,x,x^{2},\cdots,x^{n},\cdots\right\}$ 正交化得到的正交多项式就是 Chebyshev 多项式，它可表示为

$$T_{n}(x)=\cos (n \arccos x),\quad|x| \leqslant 1$$

若令 $x=\cos \theta$，则

$$T_{n}(x)=\cos n \theta,\quad 0 \leqslant \theta \leqslant \pi$$

性质

• 递推关系 $$\begin{array}{c} T_{0}(x)=1,\quad T_{1}(x)=x \\ T_{n+1}(x)=2 x T_{n}(x)-T_{n-1}(x) \quad(n=1,2,\cdots) \end{array}$$
• $T_{n}(x)$ 对零的偏差最小，它可写成如下定理
  ★ 在区间 $[-1,1]$ 上所有最高项系数为 $1$ 的一切 $n$ 次多项式中，$\omega_{n}(x) =\frac{1}{2^{n-1}} T_{n}(x)$ 与零的偏差最小，其偏差为 $\frac{1}{2^{n-1}}$
• Chebyshev 多项式 $\left\{T_{k}(x)\right\}$ 在区间 $[-1,1]$ 上带权 $\rho(x)=1 / \sqrt{1-x^{2}}$ 正交，且 $$\int_{-1}^{1} \frac{T_{n}(x) T_{m}(x) \mathrm{d} x}{\sqrt{1-x^{2}}}=\left\{\begin{array}{ll} 0,& n \neq m,\\ \frac{\pi}{2},& n=m \neq 0,\\ \pi,& n=m=0 \end{array}\right.$$ 事实上，令 $x=\cos \theta$，则 $\mathrm{d} x=-\sin \theta \mathrm{d} \theta$，于是 $$\int_{-1}^{1} \frac{T_{n}(x) T_{m}(x)}{\sqrt{1-x^{2}}} \mathrm{~d} x=\int_{0}^{\pi} \cos n \theta \cos m \theta \mathrm{~d} \theta=\left\{\begin{array}{ll} 0,& n \neq m,\\ \frac{\pi}{2},& n=m \neq 0,\\ \pi,& n=m=0 \end{array}\right.$$
• $T_{2 k}(x)$ 只含 $x$ 的偶次幂，$T_{2 k+1}(x)$ 只含 $x$ 的奇次幂。
• $T_{n}(x)$ 在区间 $[-1,1]$ 上有 $n$ 个零点 $x_{k}=\cos \frac{2 k-1}{2 n} \pi,k=1,\cdots,n$。
  此外，实际计算中时常要求 $x^{n}$ 用 $T_{0},T_{1},\cdots,T_{n}$ 的线性组合表示，其公式为 $$x^{n}=2^{1-n} \sum_{k=0}^{\left[\frac{n}{2}\right]}\binom{n}{k} T_{n-2 k}(x)$$

常用低次多项式

$$\begin{array}{l} T_{0}(x)=1 \\ T_{1}(x)=x \\ T_{2}(x)=2 x^{2}-1 \\ T_{3}(x)=4 x^{3}-3 x \\ T_{4}(x)=8 x^{4}-8 x^{2}+1 \\ T_{5}(x)=16 x^{5}-20 x^{3}+5 x \\ T_{6}(x)=32 x^{6}-48 x^{4}+18 x^{2}-1 \\ T_{7}(x)=64 x^{7}-112 x^{5}+56 x^{3}-7 x \\ T_{8}(x)=128 x^{8}-256 x^{6}+160 x^{4}-32 x^{2}+1 \end{array}$$

★求函数最佳逼近的方法

最佳零次一致逼近多项式

令 $m=\inf_{a\le x \le b} f(x)$，$M=\sup_{a\le x \le b} f(x)$，则最佳零次逼近多项式为

$$P_{0}^{*}=\frac{M-m}{2}$$

最佳一次一致逼近多项式

假定 $f(x) \in C^{2}[a,b]$，且 $f^{\prime \prime}(x)$ 在 $(a,b)$ 内不变号。

先求

$$a_{1}=\frac{f(b)-f(a)}{b-a}$$

由于 $f^{\prime \prime}(x)$ 在 $[a,b]$ 上不变号，故 $f^{\prime}(x)$ 单调，$f^{\prime}(x)-a_{1}$ 在 $(a,b)$ 内只有一个零点，解得 $x_{2}$，然后求 $a_{0}$：

$$a_{0}=\frac{f(a)+f\left(x_{2}\right)}{2}-\frac{f(b)-f(a)}{b-a} \frac{a+x_{2}}{2}\\$$

最佳一次逼近多项式 $P_{1}^{*}(x)=a_{0}+a_{1} x$，代入得

$$y=\frac{1}{2}\left[f(a)+f\left(x_{2}\right)\right]+a_{1}\left(x-\frac{a+x_{2}}{2}\right)$$

最佳 n 次平方逼近多项式

1. 写出矩阵形式法方程 $$\boldsymbol{H}\boldsymbol{a}=\boldsymbol{d}$$ 其中 $$\begin{array}{c} \boldsymbol{H}=\left(\begin{array}{cccc} 1 & 1 / 2 & \cdots & 1 /(n+1) \\ 1 / 2 & 1 / 3 & \cdots & 1 /(n+2) \\ \vdots & \vdots & & \vdots \\ 1 /(n+1) & 1 /(n+2) & \cdots & 1 /(2 n+1) \end{array}\right), \\ \boldsymbol{a}=\left(a_{0}, a_{1}, \cdots, a_{n}\right)^{\mathrm{T}}, \\ \boldsymbol{d}=\left(d_{0}, d_{1} \cdots, d_{n}\right)^{\mathrm{T}} \\ \end{array}$$ 其中 $$d_{k}=\left(f, x^{k}\right)=\int_{0}^{1} f(x) x^{k} \mathrm{~d} x\,(k=0,1, \cdots, n)$$
2. 解方程组得到 $$a_{0}^{*}, a_{1}^{*}, \cdots, a_{n}^{*}$$
3. 得到要求的多项式为 $$S^{*}(x)=a_{0}^{*}+a_{1}^{*} x+\cdots+a_{n}^{*} x^{n}$$
4. 计算平方误差 $$\|\delta\|_{2}^{2}=\|f\|_{2}^{2}-\sum_{k=0}^{n} a_{k}^{*}\left(x^k, f\right)$$
5. 计算最大误差 $$\|\delta\|_{\infty}=\max _{a \leqslant x \leqslant b}\left|f(x)-S^{*}(x)\right|$$

多项式函数的最佳 n 次一致逼近多项式

若目标区间是 $[-1,1]$，设 $f(x)$ 的最高次项为 $a_{n+1}x^{n+1}$，由切比雪夫多项式的性质，最佳 n 次逼近多项式 $P_{n}^{*}(x)$ 满足

$$f(x)-P_{n}^{*}=\frac{a_{n+1}}{2^{n}}T_{n+1}(x)$$

从而求得 $P_{n}^{*}$

若目标区间不是 $[-1,1]$（假设是 $[a,b]$），则先换元（平移+缩放）到 $[-1,1]$。

1. 令 $z=\frac{2}{b-a}x-\frac{b+a}{b-a}$，则 $z\in [-1,1]$
2. 代入 $f(x)$ 求得 $g(z)$
3. 计算 $g(z)$ 的最佳 n 次逼近多项式 $P_{n}^{*}(z)$
4. 代入 $x=\frac{b-a}{2}z+\frac{b+a}{2}$ 得到 $P_{n}^{*}(x)$

最小二乘法

定义

对于给定的一组数据 $\left(x_{i},y_{i}=f(x_{i})\right) (i=0,1,\cdots,m)$，要求在函数空间 $\Phi=\operatorname{span}\left\{\varphi_{0},\varphi_{1},\cdots,\varphi_{n}\right\}$ 中找一个函数 $y=S^{*}(x)$，使误差平方和达到最小：

$$\|\boldsymbol{\delta}\|_{2}^{2}=\sum_{i=0}^{m} \delta_{i}^{2}=\sum_{i=0}^{m}\left[S^{*}\left(x_{i}\right)-y_{i}\right]^{2}=\min _{S(x) \in \Phi} \sum_{i=1}^{m}\left[S\left(x_{i}\right)-y_{i}\right]^{2},$$

其中

$$S(x)=a_{0} \varphi_{0}(x)+a_{1} \varphi_{1}(x)+\cdots+a_{n} \varphi_{n}(x) \quad(n<m)$$

为了使问题的提法更有一般性，通常把最小二乘法中的项都考虑为加权平方和

$$\|\boldsymbol{\delta}\|_{2}^{2}=\sum_{i=0}^{m} \omega\left(x_{i}\right)\left[S\left(x_{i}\right)-f\left(x_{i}\right)\right]^{2}$$

这里 $\omega(x) \geqslant 0$ 是 $[a,b]$ 上的权函数，它表示不同点 $\left(x_{i},f\left(x_{i}\right)\right)$ 处的数据比重不同。

问题可以转化为多元函数求最值问题，求偏导得到条件：

$$\frac{\partial I}{\partial a_{k}}=2 \sum_{i=0}^{m} \omega\left(x_{i}\right)\left[\sum_{j=0}^{n} a_{j} \varphi_{j}\left(x_{i}\right)-f\left(x_{i}\right)\right] \varphi_{k}\left(x_{i}\right)=0 \quad(k=0,1,\cdots,n)$$

若记（注意和函数内积的区别！）

$$\begin{array}{c} \left(\varphi_{j},\varphi_{k}\right)=\sum_{i=0}^{m} \omega\left(x_{i}\right) \varphi_{j}\left(x_{i}\right) \varphi_{k}\left(x_{i}\right),\\ \left(f,\varphi_{k}\right)=\sum_{i=0}^{m} \omega\left(x_{i}\right) f\left(x_{i}\right) \varphi_{k}\left(x_{i}\right) \equiv d_{k} \quad(k=0,1,\cdots,n), \end{array}$$

则方程组可改写为

$$\sum_{j=0}^{n}\left(\varphi_{k},\varphi_{j}\right) a_{j}=d_{k} \quad(k=0,1,\cdots,n)$$

它也可写成矩阵形式

$$\begin{array}{c} \boldsymbol{G} \boldsymbol{a}=\boldsymbol{d},\\ \boldsymbol{a}=\left(a_{0},a_{1},\cdots,a_{n}\right)^{\mathrm{T}},\quad \boldsymbol{d}=\left(d_{0},d_{1},\cdots,d_{n}\right)^{\mathrm{T}},\\ \boldsymbol{G}=\left(\begin{array}{cccc} \left(\varphi_{0},\varphi_{0}\right) & \left(\varphi_{0},\varphi_{1}\right) & \cdots & \left(\varphi_{0},\varphi_{n}\right) \\ \left(\varphi_{1},\varphi_{0}\right) & \left(\varphi_{1},\varphi_{1}\right) & \cdots & \left(\varphi_{1},\varphi_{n}\right) \\ \vdots & \vdots & & \vdots \\ \left(\varphi_{n},\varphi_{0}\right) & \left(\varphi_{n},\varphi_{1}\right) & \cdots & \left(\varphi_{n},\varphi_{n}\right) \end{array}\right), \end{array}$$

由于 $\varphi_{0},\varphi_{1},\cdots,\varphi_{n}$ 线性无关，故 $|\boldsymbol{G}| \neq 0$，方程组存在唯一的解

$$a_{k}=a_{k}^{*} \quad(k=0,1,\cdots,n),$$

从而得到函数 $f(x)$ 的最小二乘解为

$$S^{*}(x)=a_{0}^{*} \varphi_{0}(x)+a_{1}^{*} \varphi_{1}(x)+\cdots+a_{n}^{*} \varphi_{n}(x)$$

第四章 数值积分与数值微分

常用低级积分方法

梯形公式

$$T =\frac{b-a}{2}[f(a)+f(b)]$$

矩形公式

$$R =(b-a) f\left(\frac{a+b}{2}\right)$$

代数精度

如果某个求积公式对于次数不大于 $m$ 的多项式均能准确地成立，但对于 $m+1$ 次多项式就不一定准确，则称该求积公式具有 m 次代数精度。

插值型积分

设给定一组节点

$$a \leqslant x_{0}<x_{1}<x_{2}<\cdots<x_{n} \leqslant b$$

且已知函数 $f(x)$ 在这些节点上的值，作插值函数 $L_{n}(x)$。由于代数多项式 $L_{n}(x)$ 的原函数是容易求出的，取 $I_{n}=\int_{a}^{b} L_{n}(x) \mathrm{d} x$ 作为积分 $I=\int_{a}^{b} f(x) \mathrm{d} x$ 的近似值，这样构造出的求积公式

$$I_{n}=\sum_{k=0}^{n} A_{k} f\left(x_{k}\right)$$

是插值型的，其中求积系数 $A_{k}$ 通过插值基函数 $l_{k}(x)$ 积分

$$A_{k}=\int_{a}^{b} l_{k}(x) \mathrm{d} x$$

得出

由插值余项定理即知，对于插值型的求积公式，其余项

$$R[f]=I-I_{n}=\int_{a}^{b} \frac{f^{n+1}(\xi)}{(n+1)!} w(x) \mathrm{d} x,$$

其中 $\xi$ 与变量 $x$ 有关，$w(x)=\left(x-x_{0}\right)\left(x-x_{1}\right) \cdots\left(x-x_{n}\right)$

求积公式至少有 $n$ 次代数精度的充分必要条件是，它是插值型的。

Newton-Cotes 公式

设将积分区间 $[a,b]$ 划分为 $n$ 等份，步长 $h=\frac{b-a}{n}$，选取等距节点 $x_{k}=a+k h$ 构造出的插值型求积公式

$$I_{n}=(b-a) \sum_{k=0}^{n} C_{k}^{(n)} f\left(x_{k}\right)$$

称为 Newton-Cotes 公式，其中 $C_{k}^{(n)}$ 称为 Cotes 系数

• $n=1$ 时即为梯形公式
  • 积分余项： $$R_{\mathrm{T}}=\frac{f^{\prime \prime}(\eta)}{2} \int_{a}^{b}(x-a)(x-b) \mathrm{d} x=-\frac{f^{\prime \prime}(\eta)}{12}(b-a)^{3},\quad \eta \in(a,b)$$
• $n=2$ 时是 Simpson 公式： $$S=\frac{b-a}{6}\left[f(a)+4 f\left(\frac{a+b}{2}\right)+f(b)\right]$$
  • 积分余项： $$R_{\mathrm{S}}=\frac{f^{(4)}(\eta)}{4!} \int_{a}^{b}(x-a)(x-c)^{2}(x-b) \mathrm{d} x=-\frac{b-a}{180}\left(\frac{b-a}{2}\right)^{4} f^{(4)}(\eta)$$

• $n=4$ 时称为 Cotes 公式，其形式为

  $$C=\frac{b-a}{90}\left[7 f\left(x_{0}\right)+32 f\left(x_{1}\right)+12 f\left(x_{2}\right)+32 f\left(x_{3}\right)+7 f\left(x_{4}\right)\right]$$

  这里 $x_{k}=a+k h$，$h=\frac{b-a}{4}$

  • 积分余项： $$R_{\mathrm{C}}=I-C=-\frac{2(b-a)}{945}\left(\frac{b-a}{4}\right)^{6} f^{(6)}(\eta)$$

• 当阶 $n$ 为偶数时，Newton-Cotes 公式至少有 $n+1$ 次代数精度

Cotes 系数

引进变换

$$x=a+t h,$$

则有

$$C_{k}^{(n)}=\frac{h}{b-a} \int_{0}^{n} \prod_{\substack{j=0 \\ j \neq k}}^{n} \frac{t-j}{k-j} \mathrm{~d} t=\frac{(-1)^{n-k}}{n k!(n-k)!} \int_{0}^{n} \prod_{\substack{j = 0 \\ j \neq k}}^{n}(t-j) \mathrm{d} t$$

常用 Cotes 系数：

$$\begin{array}{|c|c|c|c|c|c|c|c|c|c|} \hline n & C_{k}^{(n)} \\ \hline 1 & \frac{1}{2} & \frac{1}{2} & & & & & & & \\ \hline 2 & \frac{1}{6} & \frac{2}{3} & \frac{1}{6} & & & & & & \\ \hline 3 & \frac{1}{8} & \frac{3}{8} & \frac{3}{8} & \frac{1}{8} & & & & & \\ \hline 4 & \frac{7}{90} & \frac{16}{45} & \frac{2}{15} & \frac{16}{45} & \frac{7}{90} & & & & \\ \hline 5 & \frac{19}{288} & \frac{25}{96} & \frac{25}{144} & \frac{25}{144} & \frac{25}{96} & \frac{19}{288} & & & \\ \hline 6 & \frac{41}{840} & \frac{9}{35} & \frac{9}{280} & \frac{34}{105} & \frac{9}{280} & \frac{9}{35} & \frac{41}{840} & & \\ \hline 7 & \frac{751}{17280} & \frac{3577}{17280} & \frac{1323}{17280} & \frac{2989}{17280} & \frac{2989}{17280} & \frac{1323}{17280} & \frac{3577}{17280} & \frac{751}{17280} & \\ \hline 8 & \frac{989}{28350} & \frac{5888}{28350} & \frac{-928}{28350} & \frac{10496}{28350} & \frac{-4540}{28350} & \frac{10496}{28350} & \frac{-928}{28350} & \frac{5888}{28350} & \frac{989}{28350} \\ \hline \end{array}$$

可以看到，当 $n \geqslant 8$ 时，Cotes 系数有正有负，这时稳定性得不到保证。因此，实际计算中不用高阶的 Newton-Cotes 公式。

复化求积法

设将积分区间 $[a,b]$ 划分为 $n$ 等份，步长 $h=\frac{b-a}{n}$，分点为 $x_{k}=a+k h,k=0,1,\cdots,n$。所谓复化求积法，就是先用低阶的 Newton-Cotes 公式求得每个子区间 $\left[x_{k},x_{k+1}\right]$ 上的积分值 $I_{k}$，然后再求和，用 $\sum_{k=0}^{n-1} I_{k}$ 作为所求积分 $I$ 的近似值。

• 复化梯形公式 $$T_{n}=\sum_{k=0}^{n-1} \frac{h}{2}\left[f\left(x_{k}\right)+f\left(x_{k+1}\right)\right]=\frac{h}{2}\left[f(a)+2 \sum_{k=1}^{n-1} f\left(x_{k}\right)+f(b)\right]$$
  • 积分余项 $$I-T_{n}=\sum_{k=0}^{n-1}\left[-\frac{h^{3}}{12} f^{\prime \prime}\left(\eta_{k}\right)\right]=-\frac{b-a}{12} h^{2} f^{\prime \prime}(\eta)$$
• 复化 Simpson 公式 $$\begin{aligned} S_{n} & =\sum_{k=0}^{n-1} \frac{h}{6}\left[f\left(x_{k}\right)+4 f\left(x_{k+\frac{1}{2}}\right)+f\left(x_{k+1}\right)\right] \\ & =\frac{h}{6}\left[f(a)+4 \sum_{k=0}^{n-1} f\left(x_{k+\frac{1}{2}}\right)+2 \sum_{k=1}^{n-1} f\left(x_{k}\right)+f(b)\right] \end{aligned}$$ 其中 $x_{k+\frac{1}{2}} = \frac{1}{2}(x_{k}+x_{k+1})$
  • 积分余项 $$I-S_{n}=-\frac{b-a}{180}\left(\frac{h}{2}\right)^{4} f^{(4)}(\eta),\quad \eta \in[a,b]$$
• 复化 Cotes 公式 $$\begin{aligned} C_{n}= & \frac{h}{90}\left[7 f(a)+32 \sum_{k=0}^{n-1} f\left(x_{k+\frac{1}{4}}\right)+12 \sum_{k=0}^{n-1} f\left(x_{k+\frac{1}{2}}\right)\right. \\ & \left.+32 \sum_{k=0}^{n-1} f\left(x_{k+\frac{3}{4}}\right)+14 \sum_{k=1}^{n-1} f\left(x_{k}\right)+7 f(b)\right] \end{aligned}$$ 其中 $$\begin{array}{c} x_{k+\frac{1}{4}} = \frac{1}{4} x_{k} + \frac{3}{4} x_{k+1} \\ x_{k+\frac{1}{2}} = \frac{1}{2} x_{k} + \frac{1}{2} x_{k+1} \\ x_{k+\frac{3}{4}} = \frac{3}{4} x_{k} + \frac{1}{4} x_{k+1} \\ \end{array}$$
  • 积分余项 $$I-C_{n}=-\frac{2(b-a)}{945}\left(\frac{h}{4}\right)^{6} f^{(6)}(\eta),\quad \eta \in[a,b]$$

Romberg 算法

变步长二分递推

设将求积区间 $[a,b]$ 分成 $n$ 等份，则一共有 $n+1$ 个分点，按梯形公式计算积分值 $T_{n}$，需要提供 $n+1$ 个函数值。

如果将求积区间再二分一次，则分点增至 $2 n+1$ 个。

将二分前后两个积分值联系起来加以考察，注意到每个子区间 $\left[x_{k},x_{k+1}\right]$ 经过二分只增加了一个分点 $x_{k+\frac{1}{2}}=\frac{1}{2}\left(x_{k}+x_{k+1}\right)$，用复化梯形公式求得该子区间上的积分值为

$$\frac{h}{4}\left[f\left(x_{k}\right)+2 f\left(x_{k+\frac{1}{2}}\right)+f\left(x_{k+1}\right)\right]$$

注意，这里 $h=\frac{b-a}{n}$ 代表二分前的步长。将每个子区间上的积分值相加，得

$$T_{2 n}=\frac{h}{4} \sum_{k=0}^{n-1}\left[f\left(x_{k}\right)+f\left(x_{k+1}\right)\right]+\frac{h}{2} \sum_{k=0}^{n-1} f\left(x_{k+\frac{1}{2}}\right),$$

从而可导出下列递推公式：

$$T_{2 n}=\frac{1}{2} T_{n}+\frac{h}{2} \sum_{k=0}^{n-1} f\left(x_{k+\frac{1}{2}}\right)$$

估计补偿

根据梯形法的误差公式，积分值 $T_{n}$ 的截断误差大致与 $h^{2}$ 成正比，因此当步长二分后，截断误差将减至原有误差的 $1 / 4$，即有

$$\frac{I-T_{2 n}}{I-T_{n}} \approx \frac{1}{4}$$

将上式移项整理，可得

$$I-T_{2 n} \approx \frac{1}{3}\left(T_{2 n}-T_{n}\right)$$

由此可见，积分近似值 $T_{2 n}$ 的误差大致等于 $\frac{1}{3}\left(T_{2 n}-T_{n}\right)$，因此，如果用这个误差值作为 $T_{2 n}$ 的一种补偿，可以期望更好的结果：

$$\bar{T}=T_{2 n}+\frac{1}{3}\left(T_{2 n}-T_{n}\right)=\frac{4}{3} T_{2 n}-\frac{1}{3} T_{n}$$

而事实上这与 Simpson 法的积分值相等

$$S_n = \frac{4}{3} T_{2 n}-\frac{1}{3} T_{n} = \bar{T}$$

同样，$S_{n}$ 补偿之后变成了 Cotes 公式

$$C_{n} =\frac{16}{15} S_{2 n}-\frac{1}{15} S_{n}$$

$C_{n}$ 补偿之后变成了 Romberg 公式

$$R_{n} =\frac{64}{63} C_{2 n}-\frac{1}{63} C_{n}$$

Richardson 外推加速法

（理论依据）设 $f(x) \in C^{\infty}[a,b]$，则成立

$$T(h)=I+a_{1} h^{2}+a_{2} h^{4}+a_{3} h^{6}+\cdots+a_{k} h^{2 k}+\cdots,$$

其中系数 $a_{k}(k=1,2,\cdots)$ 与 $h$ 无关。

设以 $T_{0}^{(k)}$ 表示二分 $k$ 次后求得的梯形值，且以 $T_{m}^{(k)}$ 表示序列 $\left\{T_{0}^{(k)}\right\}$ 的 $m$ 次加速值，$h_{}$则依递推公式，即

$$T_{m}^{(k)}=\frac{4^{m}}{4^{m}-1} T_{m-1}^{(k+1)}-\frac{1}{4^{m}-1} T_{m-1}^{(k)} \quad(k=1,2,\cdots)$$

可以逐行构造出T 数表：

$$\begin{array}{ccccc} T_{0}^{(0)} & & & & \\ T_{0}^{(1)} & T_{1}^{(0)} & & & \\ T_{0}^{(2)} & T_{1}^{(1)} & T_{2}^{(0)} & & \\ T_{0}^{(3)} & T_{1}^{(2)} & T_{2}^{(1)} & T_{3}^{(0)} & \\ \vdots & \vdots & \vdots & \vdots & \ddots \end{array}$$

可以证明，如果 $f(x)$ 充分光滑，那么 $T$ 数表每一列的元素及对角线元素均收敛到所求的积分值 $I$，即

$$\lim _{k \rightarrow \infty} T_{m}^{(k)}=I \quad(m \mathrm{固定}),\quad \lim _{m \rightarrow \infty} T_{m}^{(0)}=I$$

• 初值： $$T_0^{(0)} = \frac{b-a}{2}\left(f(a)+f(b)\right)$$
• 纵向递推： $$T_{0}^{(k)}=\frac{1}{2} T_{0}^{(k-1)}+\frac{b-a}{2^{k}} \sum_{i=0}^{2^{k-1}-1} f\left(x_{i+\frac{1}{2}}\right)$$
• 横向递推： $$T_{m}^{(k)}=\frac{4^{m}}{4^{m}-1} T_{m-1}^{(k+1)}-\frac{1}{4^{m}-1} T_{m-1}^{(k)}$$

Gauss 公式

定义

在机械求积公式

$$\int_{a}^{b} f(x) \mathrm{d} x \approx \sum_{k=0}^{n} A_{k} f\left(x_{k}\right)$$

中含有 $2 n+2$ 个待定参数 $x_{k},A_{k}(k=0,1,\cdots,n)$，适当选择这些参数，有可能使求积公式具有 $2 n+1$ 次代数精度，这类求积公式称为 Gauss 公式

Gauss 点

如果求积公式具有 $2 n+1$ 次代数精度，则称其节点 $x_{k}(k=0,1,\cdots,n)$ 是 Gauss 点

• 对于插值型求积公式，其节点 $x_{k}(k=0,1,\cdots,n)$ 是 Gauss 点的充分必要条件，是以这些点为零点的多项式 $\omega(x)=\prod_{k=0}^{n}\left(x-x_{k}\right)$ 与任意次数不超过 $n$ 的多项式 $P(x)$ 均正交，即 $$\int_{a}^{b} P(x) \omega(x) \mathrm{d} x=0$$

Gauss-Legendre 公式

不失一般性，可取 $a=-1,b=1$ 而考察区间 $[-1,1]$ 上的 Gauss 公式。显然，Legendre 多项式 $P_{n+1}(x)$ 的零点就是求积公式的 Gauss 点。这样子的 Gauss 公式特别地称为 Gauss-Legendre 公式。

常用低次 Gauss-Legendre 公式

• $n=0$：

  $$\int_{-1}^{1} f(x) \mathrm{d} x \approx 2 f(0)$$

• $n=1$：

  $$\int_{-1}^{1} f(x) \mathrm{d} x \approx f(-\frac{1}{\sqrt{3}}) + f(\frac{1}{\sqrt{3}})$$

• $n=2$：

  $$\int_{-1}^{1} f(x) \mathrm{d} x \approx \frac{5}{9} f(-\frac{\sqrt{15}}{5}) + \frac{8}{9} f(0) + \frac{5}{9} f(\frac{\sqrt{15}}{5})$$

• $n=3$：

  $$\int_{-1}^{1} f(x) \mathrm{d} x \approx \frac{18+\sqrt{30}}{36} \Sigma f(\pm \sqrt{\frac{3}{7} - \frac{2}{7}\sqrt{\frac{6}{5}}})+ \frac{18-\sqrt{30}}{36} \Sigma f(\pm \sqrt{\frac{3}{7} + \frac{2}{7}\sqrt{\frac{6}{5}}})$$

• $n=4$：

  $$\int_{-1}^{1} f(x) \mathrm{d} x \approx \frac{128}{225} f(0)+ \frac{322+13\sqrt{70}}{900} \Sigma f(\pm \frac{1}{3}\sqrt{5 - 2\sqrt{\frac{10}{7}}})+ \frac{322-13\sqrt{70}}{900} \Sigma f(\pm \frac{1}{3}\sqrt{5 + 2\sqrt{\frac{10}{7}}})$$

收敛性

对于 Gauss 公式，其余项

$$R(x)=\int_{a}^{b} f(x) \mathrm{d} x-\sum_{k=0}^{n} A_{k} f\left(x_{k}\right)=\frac{f^{(2 n+2)}(\xi)}{(2 n+2)!} \int_{a}^{b} \omega^{2}(x) \mathrm{d} x,$$

这里 $\omega(x)=\left(x-x_{0}\right)\left(x-x_{1}\right) \cdots\left(x-x_{n}\right)$

稳定性

Gauss 公式的求积系数 $A_{k}(k=0,1,\cdots,n)$ 全是正的

复化两点 Gauss 型积分公式

将区间 $[a,b]$ 进行 $m$ 等分，每个小区间 $x \in\left[x_{k},x_{k+1}\right]$ 上的两点 Gauss 积分为

$$\frac{h}{2}\left(f\left(x_{k,0}\right)+f\left(x_{k,1}\right)\right)$$

其中

$$\begin{array}{c} h=\frac{b-a}{m},\\ x_{k,0}=\frac{x_{k}+x_{k+1}}{2}-\frac{h}{2} \sqrt{\frac{1}{3}},\\ x_{k,1}=\frac{x_{k}+x_{k+1}}{2}+\frac{h}{2} \sqrt{\frac{1}{3}} \\ \end{array}$$

求和得到

$$G_{m}=\sum_{k=0}^{m-1} \frac{h}{2}\left(f\left(x_{k,0}\right)+f\left(x_{k,1}\right)\right)$$

第五章 常微分方程数值解法

问题背景

一阶方程的初值问题：

$$\left\{\begin{array}{l} y^{\prime}=f(x,y),\\ y\left(x_{0}\right)=y_{0} \end{array}\right.$$

注意：

• 在第一个方程中，$x$ 和 $y$ 表示所有可能的点，可以理解成参数方程 $x(t)$ 和 $y(t)$，而 $f$ 则是关于 $x(t)$ 和 $y(t)$ 的函数
• 在第二个方程中 $y(x_0)$ 表示与 $x_0$ 对应的 $y$ 值，这个方程定义了方程的初始条件
• 在 $O x y$ 平面上，微分方程的解 $y=y(x)$ 称为它的积分曲线

所谓数值解法，就是寻求解 $y(x)$ 在一系列离散节点

$$x_{1}<x_{2}<\cdots<x_{n}<x_{n+1}<\cdots$$

上的近似值 $y_{1},y_{2},\cdots,y_{n},y_{n+1},\cdots$，相邻两个节点的间距 $h=x_{n+1}-x_{n}$ 称为步长，$x_{n}=x_{0}+n h,n=0,1,2,\cdots$

初值问题的数值解法有个基本特点，它们都采取“步进式”，即求解过程顺着节点排列的次序一步一步地向前推进。描述这类算法，只要给出用已知信息 $y_{n},y_{n-1},y_{n-2},\cdots$ 计算 $y_{n+1}$ 的递推公式即可。这种计算公式称为差分格式

前进欧拉

由方程可知，积分曲线上一点 $(x_n,y_n)$ 的切线斜率等于函数 $f(x,y)$ 在 $(x_n,y_n)$ 的值，故可以将此处的斜率作为小区间 $[x_n, x_{n+1}]$ 内的斜率（即用一条直线来近似积分曲线），于是可以写出直线方程：

$$\frac{y_{n+1}-y_{n}}{x_{n+1}-x_{n}}=f\left(x_{n},y_{n}\right)$$

那么 $y_{n+1}$ 的值为

$$y_{n+1}=y_{n}+h f\left(x_{n},y_{n}\right)$$

这种方法叫做前进欧拉格式（简称欧拉格式）

截断误差：

$$y\left(x_{n+1}\right)-y_{n+1}=\frac{h^{2}}{2} y^{\prime \prime}(\xi) \approx \frac{h^{2}}{2} y^{\prime \prime}\left(x_{n}\right)$$

后退欧拉

若使用向后的差商代替方程中的导数，即

$$\frac{y_{n}-y_{n-1}}{x_{n}-x_{n-1}}=f\left(x_{n},y_{n}\right)$$

移项，$n=n+1$，得到

$$y_{n+1}=y_{n}+h f\left(x_{n+1},y_{n+1}\right)$$

这种方法叫做后退欧拉格式，这种格式实际上是关于 $y_{n+1}$ 的一个方程，因此是隐式的，需要使用迭代法求解。

迭代公式：

$$\left\{\begin{array}{l} y_{n+1}^{(0)}=y_{n}+h f\left(x_{n},y_{n}\right) \\ y_{n+1}^{(k+1)}=y_{n}+h f\left(x_{n+1},y_{n+1}^{(k)}\right) \end{array}\right.$$

截断误差：

$$y\left(x_{n+1}\right)-y_{n+1} \approx-\frac{h^{2}}{2} y^{\prime \prime}\left(x_{n}\right)$$

梯形公式

比较前进欧拉格式与后退欧拉格式的误差公式，可以看到，如果将这两种方法进行算术平均，即可消除误差的主要部分 $\pm \frac{h^{2}}{2} y^{\prime \prime}_{n}$ 而获得更高的精度。这种平均化方法通常称为梯形方法，其计算格式为

$$y_{n+1}=y_{n}+\frac{h}{2}\left[f\left(x_{n},y_{n}\right)+f\left(x_{n+1},y_{n+1}\right)\right]$$

迭代公式：

$$\left\{\begin{array}{l} y_{n+1}^{(0)}=y_{n}+h f\left(x_{n},y_{n}\right),\\ y_{n+1}^{(k+1)}=y_{n}+\frac{h}{2}\left[f\left(x_{n},y_{n}\right)+f\left(x_{n+1},y_{n+1}^{(k)}\right)\right] \end{array}(k=0,1,2,\cdots)\right.$$

改进的欧拉格式

先用欧拉格式求得一个初步的近似值 $\bar{y}_{n+1}$，称之为预测值。预测值 $\bar{y}_{n+1}$ 的精度可能很差，再用梯形公式将它校正一次，得 $y_{n+1}$，这个结果称为校正值。而这样建立的预测-校正系统通常称为改进的欧拉格式：

$$\left\{\begin{array}{c} \bar{y}_{n+1}=y_{n}+h f\left(x_{n},y_{n}\right),\\ y_{n+1}=y_{n}+\frac{h}{2}\left[f\left(x_{n},y_{n}\right)+f\left(x_{n+1},\bar{y}_{n+1}\right)\right] \end{array}\right.$$

这一计算格式亦可表示为

$$y_{n+1}=y_{n}+\frac{h}{2}\left[f\left(x_{n},y_{n}\right)+f\left(x_{n}+h,y_{n}+h f\left(x_{n},y_{n}\right)\right)\right],$$

或表示为下列平均化形式：

$$\left\{\begin{array}{l} y_{p}=y_{n}+h f\left(x_{n},y_{n}\right) \\ y_{c}=y_{n}+h f\left(x_{n+1},y_{p}\right) \\ y_{n+1}=\frac{1}{2}\left(y_{p}+y_{c}\right) \end{array}\right.$$

欧拉两步格式

$$y_{n+1}=y_{n-1}+2 h f\left(x_{n},y_{n}\right)$$

迭代方程：

$$\begin{array}{c} \bar{y}_{n+1}=y_{n-1}+2 h f\left(x_{n},y_{n}\right) \\ y_{n+1}=y_{n}+\frac{h}{2}\left[f\left(x_{n},y_{n}\right)+f\left(x_{n+1},\bar{y}_{n+1}\right)\right] \end{array}$$

Taylor 格式

设初值问题有解 $y=y(x)$，用 Taylor 展开，有

$$y\left(x_{n+1}\right)=y\left(x_{n}\right)+h y^{\prime}\left(x_{n}\right)+\frac{h^{2}}{2!} y^{\prime \prime}\left(x_{n}\right)+\frac{h^{3}}{3!} y^{\prime \prime \prime}\left(x_{n}\right)+\cdots,$$

其中 $y(x)$ 的各阶导数依据所给方程可以用函数 $f$ 来表达。下面引进函数序列 $f^{(j)}(x,y)$ 来描述求导过程，即

$$y^{\prime}=f \equiv f^{(0)},\quad y^{\prime \prime}=\frac{\partial f^{(0)}}{\partial x}+f \frac{\partial f^{(0)}}{\partial y} \equiv f^{(1)},\quad y^{\prime \prime \prime}=\frac{\partial f^{(1)}}{\partial x}+f \frac{\partial f^{(1)}}{\partial y} \equiv f^{(2)},\quad \cdots$$

一般地

$$y^{(j)}=\frac{\partial f^{(j-2)}}{\partial x}+f \frac{\partial f^{(j-2)}}{\partial y} \equiv f^{(j-1)},\quad j=2,3,\cdots,$$

在展开式的右端截取若干项，并且在 $\left(x_{n},y_{n}\right)$ 计算系数 $y^{(j)}\left(x_{n}\right)$ 的近似值 $y_{n}^{(j)}$，结果导出下列 Taylor 格式

$$y_{n+1}=y_{n}+h y_{n}^{\prime}+\frac{h^{2}}{2!} y_{n}^{\prime \prime}+\cdots+\frac{h^{p}}{p!} y_{n}^{(p)}$$

注意此处的求导可以看成参数方程 $y(t)$ 对 $t$ 求导，而不是 $y$ 对 $x$ 求导！

提高 Taylor 格式的阶 $p$ 即可提高计算结果的精度。显然，$p$ 阶 Taylor 格式的局部截断误差为

$$y\left(x_{n+1}\right)-y_{n+1}=\frac{h^{p+1}}{(p+1)!} y^{(p+1)}(\xi),\quad x_{n}<\xi<x_{n+1}$$

• 如果一种方法的局部截断误差为 $O\left(h^{p+1}\right)$，则称该方法具有 p 阶精度

Runge-Kutta 方法

考察差商 $\frac{y\left(x_{n+1}\right)-y\left(x_{n}\right)}{h}$，根据微分中值定理，存在 $0<\theta<1$，使得

$$\frac{y\left(x_{n+1}\right)-y\left(x_{n}\right)}{h}=y^{\prime}\left(x_{n}+\theta h\right),$$

于是，利用所给方程 $y^{\prime}=f(x,y)$ 得到

$$y\left(x_{n+1}\right)=y\left(x_{n}\right)+h f\left(x_{n}+\theta h,y\left(x_{n}+\theta h\right)\right)$$

设 $K^{*}=f\left(x_{n}+\theta h,y\left(x_{n}+\theta h\right)\right)$，称 $K^{*}$ 为区间 $\left[x_{n},x_{n+1}\right]$ 上的平均斜率。由此可见，只要对平均斜率提供一种算法，那么就相应地导出一种计算格式。

欧拉格式简单地取点 $x_{n}$ 的斜率值 $K_{1}=f\left(x_{n},y_{n}\right)$ 作为平均斜率 $K^{*}$，精度自然很低。这个处理过程启示我们，如果设法在 $\left(x_{n},x_{n+1}\right)$ 内多预测几个点的斜率值，然后将它们加权平均作为平均斜率 $K^{*}$，则有可能构造出具有更高精度的计算格式。

二阶格式

首先推广改进的欧拉方法。随意考察区间 $\left(x_{n},x_{n+1}\right)$ 内一点

$$x_{n+p}=x_{n}+p h,\quad 0<p \leqslant 1,$$

希望用 $x_{n}$ 和 $x_{n+p}$ 两个点的斜率值 $K_{1}$ 和 $K_{2}$ 线性组合得到平均斜率 $K^{*}$，即令

$$y_{n+1}=y_{n}+h\left(\lambda_{1} K_{1}+\lambda_{2} K_{2}\right),$$

其中 $\lambda_{1},\lambda_{2}$ 为待定系数。同改进的欧拉格式一样，这里仍取 $K_{1}=f\left(x_{n},y_{n}\right)$。仿照改进的欧拉格式，先用欧拉格式提供 $y\left(x_{n+p}\right)$ 的预测值，即

$$y_{n+p}=y_{n}+p h K_{1}$$

然后再用预测值 $y_{n+p}$ 通过计算 $f$ 产生斜率值

$$K_{2}=f\left(x_{n+p},y_{n+p}\right)$$

这样设计出的计算格式具有如下形式：

$$\left\{\begin{array}{l} y_{n+1}=y_{n}+h\left(\lambda_{1} K_{1}+\lambda_{2} K_{2}\right),\\ K_{1}=f\left(x_{n},y_{n}\right) \\ K_{2}=f\left(x_{n+p},y_{n}+p h K_{1}\right) \end{array}\right.$$

欲使上式具有二阶精度，只要成立

$$\left\{\begin{array}{l} \lambda_{1}+\lambda_{2}=1,\\ \lambda_{2} p=\frac{1}{2} \end{array}\right.$$

满足条件格式统称为二阶 Runge-Kutta 格式

• 一种特殊的二阶 Runge-Kutta 格式是所谓变形的 Euler 格式，其形式如下： $$\left\{\begin{array}{l} y_{n+1}=y_{n}+h K_{2} \\ K_{1}=f\left(x_{n},y_{n}\right) \\ K_{2}=f\left(x_{n}+\frac{h}{2},y_{n}+\frac{h}{2} K_{1}\right) \end{array}\right.$$

三阶格式

为了进一步提高精度，设除 $x_{n},x_{n+p}$ 外再考察一点 $x_{n+p}=x_{n}+q h$，$p \leqslant q \leqslant 1$，并用三个点 $x_{n},x_{n+p},x_{n+q}$ 的斜率值 $K_{1},K_{2},K_{3}$ 线性组合得到平均斜率 $K^{*}$，这时计算格式

$$\left\{\begin{array}{l} y_{n+1}=y_{n}+h\left(\lambda_{1} K_{1}+\lambda_{2} K_{2}+\lambda_{3} K_{3}\right) \\ K_{1}=f\left(x_{n},y_{n}\right) \\ K_{2}=f\left(x_{n}+p h,y_{n}+p h K_{1}\right) \\ K_{3}=f\left(x_{n}+q h,y_{n}+q h\left(r K_{1}+s K_{2}\right)\right) \end{array}\right.$$

为了保证它能与三阶 Taylor 格式具有同等精度，只要满足下列方程组：

$$\left\{\begin{array}{l} r+s=1\\ \lambda_{1} + \lambda_{2} + \lambda_{3} = 1 \\ \lambda_{2} p+\lambda_{3} q=\frac{1}{2} \\ \lambda_{2} p^{2}+\lambda_{3} q^{2}=\frac{1}{3} \\ \lambda_{3} p q s=\frac{1}{6} \end{array}\right.$$

满足条件的一族格式统称三阶 Runge-Kutta 格式。

• 下面是一个特例： $$\left\{\begin{array}{l} y_{n+1}=y_{n}+\frac{h}{6}\left(K_{1}+4 K_{2}+K_{3}\right) \\ K_{1}=f\left(x_{n},y_{n}\right) \\ K_{2}=f\left(x_{n}+\frac{h}{2},y_{n}+\frac{h}{2} K_{1}\right) \\ K_{3}=f\left(x_{n}+h,y_{n}-h K_{1}+2 h K_{2}\right) \end{array}\right.$$

四阶格式

• 一个特例： $$\left\{\begin{array}{l} y_{n+1}=y_{n}+\frac{h}{6}\left(K_{1}+2 K_{2}+2 K_{3}+K_{4}\right) \\ K_{1}=f\left(x_{n},y_{n}\right) \\ K_{2}=f\left(x_{n}+\frac{h}{2},y_{n}+\frac{h}{2} K_{1}\right) \\ K_{3}=f\left(x_{n}+\frac{h}{2},y_{n}+\frac{h}{2} K_{2}\right) \\ K_{4}=f\left(x_{n}+h,y_{n}+h K_{3}\right) \end{array}\right.$$

方程组的求解

考察一阶方程组

$$y_{i}^{\prime}=f_{i}\left(x,y_{1},y_{2},\cdots,y_{N}\right) \quad(i=1,2,\cdots,N)$$

的初值问题，初始条件为

$$y_{i}\left(x_{0}\right)=y_{i}^{0} \quad(i=1,2,\cdots,N)$$

若采用向量的记号，记

$$\begin{array}{c} \boldsymbol{y}=\left(y_{1},y_{2},\cdots,y_{N}\right)^{\mathrm{T}},\quad \boldsymbol{y}^{0}=\left(y_{1}^{0},y_{2}^{0},\cdots,y_{N}^{0}\right)^{\mathrm{T}} \\ \boldsymbol{f}=\left(f_{1},f_{2},\cdots,f_{N}\right)^{\mathrm{T}} \end{array}$$

则上述方程组的初值问题可表示为

$$\left\{\begin{array}{l} \boldsymbol{y}^{\prime}=\boldsymbol{f}(x,\boldsymbol{y}) \\ \boldsymbol{y}\left(x_{0}\right)=\boldsymbol{y}_{0} \end{array}\right.$$

• Euler 法 $$\boldsymbol{y}_{n+1}=\boldsymbol{y}_{n}+h \boldsymbol{f}\left(x_{n},\boldsymbol{y}_{n}\right)$$
• 后退 Euler 法 $$\boldsymbol{y}_{n+1}=\boldsymbol{y}_{n}+h \boldsymbol{f}\left(x_{n+1},\boldsymbol{y}_{n+1}\right)$$
• 梯形格式 $$\boldsymbol{y}_{n+1}=\boldsymbol{y}_{n}+\frac{h}{2}\left(\boldsymbol{f}\left(x_{n},\boldsymbol{y}_{n}\right)+\boldsymbol{f}\left(x_{n+1},\boldsymbol{y}_{n+1}\right)\right)$$

高阶微分方程（或方程组）

原则上总可以归结为一阶方程组来求解，例如，考察下列 $m$ 阶微分方程：

$$y^{(m)}=f\left(x,y,y^{\prime},\cdots,y^{(m-1)}\right)$$

初始条件为

$$y\left(x_{0}\right)=y_{0},\quad y^{\prime}\left(x_{0}\right)=y_{0}^{\prime},\quad \cdots,\quad y^{(m-1)}\left(x_{0}\right)=y_{0}^{(m-1)}$$

只要引进新的变量 $y_{1}=y,y_{2}=y^{\prime},\cdots,y_{m}=y^{(m-1)}$，即可将 $m$ 阶方程化为如下的一阶方程组：

$$\left\{\begin{array}{l} y_{1}^{\prime} =y_{2} \\ y_{2}^{\prime} =y_{3} \\ \quad\vdots \\ y_{m-1}^{\prime} =y_{m} \\ y_{m}^{\prime} =f\left(x,y_{1},y_{2},\cdots,y_{m}\right) \end{array}\right.$$

初始条件则相应地化为

$$y_{1}\left(x_{0}\right)=y_{0},\quad y_{2}\left(x_{0}\right)=y_{0}^{\prime},\quad \cdots,\quad y_{m}\left(x_{0}\right)=y_{0}^{(m-1)}$$

特别地，对于下列二阶方程的初值问题：

$$\left\{\begin{array}{l} y^{\prime \prime}=f\left(x,y,y^{\prime}\right) \\ y\left(x_{0}\right)=y_{0},\quad y^{\prime}\left(x_{0}\right)=y_{0}^{\prime} \end{array}\right.$$

引进新的变量 $z=y^{\prime}$，即可化为下列一阶方程组的初值问题：

$$\left\{\begin{array}{l} y^{\prime}=z,\quad y\left(x_{0}\right)=y_{0} \\ z^{\prime}=f(x,y,z),\quad z\left(x_{0}\right)=y_{0}^{\prime} \end{array}\right.$$

边值问题

对于高阶微分方程，定解条件通常有两种给法

• 一种是给出了积分曲线在初始时刻的性态，这类条件称为初始条件，相应的定解问题称为初值问题
• 另一种是给出了积分曲线首末两端的性态，这类条件则称为边界条件，相应的定解问题称为边值问题

打靶法

以二阶方程

$$y^{\prime \prime}=f\left(x,y,y^{\prime}\right)$$

为例，设在区间 $a<x<b$ 上求解方程，则边界条件可以给为

$$y(a)=\alpha,\quad y(b)=\beta$$

反复调整初始时刻的斜率值 $y^{\prime}(a)=m$，使初值问题的积分曲线 $y=y(x)$ 能“命中”$y(b)=\beta$
迭代过程：

$$m_{3}=m_{1}+\frac{m_{2}-m_{1}}{\beta_{2}-\beta_{1}}\left(\beta-\beta_{1}\right)$$

差分方程法

如果所给方程是如下形式的线性方程：

$$y^{\prime \prime}+p(x) y^{\prime}+q(x) y=r(x)$$

设将积分区间 $[a,b]$ 划分为 $N$ 等份，步长 $h=\frac{b-a}{N}$，节点 $x_{n}=x_{0}+n h,n=0,1,\cdots,N$

用差商替代相应的导数，可将边值问题离散化得下列计算公式：

$$\left\{\begin{array}{l} \frac{y_{n+1}-2 y_{n}+y_{n-1}}{h^{2}}=f\left(x_{n},y_{n},\frac{y_{n+1}-y_{n-1}}{2 h}\right), \,(n=1,2,\cdots,N-1) \\ y_{0}=\alpha,\quad y_{N}=\beta \end{array}\right.$$

化简得到

$$\frac{y_{n+1}-2 y_{n}+y_{n-1}}{h^{2}}+p(x_{n}) \frac{y_{n+1}-y_{n-1}}{2 h}+q(x_{n}) y_{n}=r(x_{n}) \quad(n=1,\cdots,N-1)$$

利用边界条件消去式中的 $y_{0}$ 和 $y_{N}$，整理得到关于 $y_{n}(1 \leqslant n \leqslant N-1)$ 的下列方程组：

$$\left\{\begin{array}{l} \left(-2+h^{2} q(x_{1})\right) y_{1}+\left(1+\frac{h}{2} p(x_{1})\right) y_{2}=h^{2} r(x_{1})-\left(1-\frac{h}{2} p(x_{1})\right) \alpha \\ \left(1-\frac{h}{2} p(x_{n})\right) y_{n-1}+\left(-2+h^{2} q(x_{n})\right) y_{n}+\left(1+\frac{h}{2} p(x_{n})\right) y_{n+1}=h^{2} r(x_{n}) \quad(2 \leqslant n \leqslant N-2),\\ \left(1-\frac{h}{2} p(x_{N-1})\right) y_{N-2}+\left(-2+h^{2} q(x_{N-1})\right) y_{N-1}=h^{2} r(x_{N-1})-\left(1+\frac{h}{2} p(x_{N-1})\right) \beta \end{array}\right.$$

这样归结出的方程组是所谓三对角型的，即

$$\left(\begin{array}{ccccc} -2+h^{2} q(x_{1}) & 1+\frac{h}{2} p(x_{1}) & & & \\ 1-\frac{h}{2} p(x_{2}) & -2+h^{2} q(x_{2}) & 1+\frac{h}{2} p(x_{2}) & & \\ & \ddots & \ddots & \ddots & \\ & & 1-\frac{h}{2} p(x_{N-2}) & -2+h^{2} q(x_{N-2}) & 1+\frac{h}{2} p(x_{N-2}) \\ & & & 1-\frac{h}{2} p(x_{N-1}) & -2+h^{2} q(x_{N-1}) \end{array}\right) \boldsymbol{y}= \left(\begin{array}{c} h^{2} r(x_{1})-\left(1-\frac{h}{2} p(x_{1})\right) \alpha \\ h^{2} r(x_{2}) \\ \vdots \\ h^{2} r(x_{N-2}) \\ h^{2} r(x_{N-1})-\left(1+\frac{h}{2} p(x_{N-1})\right) \beta \end{array}\right)$$

使用追赶法求解即可。

有时边界条件按以下方式给出：

• 当 $x=a$ 时, $y^{\prime}=\alpha_{0} y+\beta_{0}$
• 当 $x=b$ 时, $y^{\prime}=\alpha_{1} y+\beta_{1}$

这里 $\alpha_{0},\beta_{0},\alpha_{1},\beta_{1}$ 均为已知常数

这时，边界条件中所包含的微商也要替换成相应的差商

$$\frac{y_{1}-y_{0}}{h}=\alpha_{0} y_{0}+\beta_{0},\quad \frac{y_{N}-y_{N-1}}{h}=\alpha_{1} y_{N}+\beta_{1}$$

它们和差分方程一起，仍然构成包含 $N+1$ 个未知数的线性方程组

• 差分问题式的解存在并且是唯一的。
• 设 $y_{n}$ 是差分问题式的（数值）解，而 $y\left(x_{n}\right)$ 是边值问题的（精确）解 $y(x)$ 在节点 $x_{n}$ 的值，则截断误差 $e_{n}=y\left(x_{n}\right)-y_{n}$ 有下列估计式： $$\left|e_{n}\right| \leqslant \frac{M(b-a)^{2}}{96} h^{2},\quad M=\max _{a \leqslant x \leqslant b}\left|y^{(4)}(x)\right|$$

第九章 矩阵的特征值与特征向量计算

基础概念

略

基本定理

如果 $\lambda_{i}(i=1,2,\cdots,n)$ 是矩阵 $\boldsymbol{A}$ 的特征值，则有

• $$\sum_{i=1}^{n} \lambda_{i}=\sum_{i=1}^{n} a_{i i}=\operatorname{tr} \boldsymbol{A}$$
• $$\operatorname{det} \boldsymbol{A}=\lambda_{1} \lambda_{2} \cdots \lambda_{n}$$

设 $\boldsymbol{A}$ 与 $\boldsymbol{B}$ 为相似矩阵（即存在非奇异阵 $\boldsymbol{T}$ 使 $\boldsymbol{B}=\boldsymbol{T}^{-1} \boldsymbol{A} \boldsymbol{T}$），则

• $\boldsymbol{A}$ 与 $\boldsymbol{B}$ 有相同的特征值；
• 若 $\boldsymbol{x}$ 是 $\boldsymbol{B}$ 的一个特征向量，则 $\boldsymbol{T} \boldsymbol{x}$ 是 $\boldsymbol{A}$ 的特征向量

（Gerschgorin’s 定理）设 $\boldsymbol{A}=\left(a_{i j}\right)_{n \times n}$，则 $\boldsymbol{A}$ 的每一个特征值必属于下述某个圆盘之中：

$$\left|\lambda-a_{i i}\right| \leqslant \sum_{\substack{j=1 \\ j \neq i}}^{n}\left|a_{i j}\right| \quad(i=1,2,\cdots,n)$$

设 $\boldsymbol{A}$ 为 $n$ 阶实对称矩阵，对于任一非零向量 $\boldsymbol{x}$，称 $R(\boldsymbol{x})=\frac{(\boldsymbol{A} \boldsymbol{x},\boldsymbol{x})}{(\boldsymbol{x},\boldsymbol{x})}$ 为对应于向量 $\boldsymbol{x}$ 的 Rayleigh 商

设 $A \in \mathbf{R}^{n \times n}$ 为对称矩阵（其特征值次序记作 $\lambda_{1} \geqslant \lambda_{2} \geqslant \cdots \geqslant \lambda_{n}$，对应的特征向量 $\boldsymbol{x}_{1},\boldsymbol{x}_{2},\cdots,\boldsymbol{x}_{n}$ 组成规范化正交组，即 $\left(\boldsymbol{x},\boldsymbol{x}_{j}\right)=\delta_{i j}$，则

• $$\lambda_{n} \leqslant \frac{(\boldsymbol{A} \boldsymbol{x},\boldsymbol{x})}{(\boldsymbol{x},\boldsymbol{x})} \leqslant \lambda_{1} \quad \text{（对于任何非零} x \in \mathbf{R}^{n} \text{）}$$
• $$\lambda_{1}=\max _{\substack{x \in \mathrm{R}^{n} \\ x \neq 0}} \frac{(\boldsymbol{A} \boldsymbol{x},\boldsymbol{x})}{(\boldsymbol{x},\boldsymbol{x})}$$
• $$\lambda_{n}=\min _{\substack{x \in \mathbf{R}^{n} \\ x \neq 0}} \frac{(\boldsymbol{A} x,\boldsymbol{x})}{(\boldsymbol{x},\boldsymbol{x})}$$

幂法

• 在一些工程、物理问题中，通常只需要求出矩阵的按模最大的特征值（称为 $\boldsymbol{A}$ 的主特征值）和相应的特征向量
• 幂法是一种计算实矩阵 $\boldsymbol{A}$ 的主特征值的一种迭代法，它最大的优点是方法简单，对于稀疏矩阵较合适，但有时收敛速度很慢

设 $\boldsymbol{A} \in \mathbf{R}^{n \times n}$ 有 $n$ 个线性无关的特征向量，主特征值 $\lambda_{1}$ 满足 $\left|\lambda_{1}\right|>\left|\lambda_{2}\right| \geqslant\left|\lambda_{3}\right| \geqslant \cdots \geqslant\left|\lambda_{n}\right|$，则对于任意非零初始向量 $\boldsymbol{v}_{0}=\boldsymbol{u}_{0}\left(a_{1} \neq 0\right)$，按下述方法构造的向量序列

$$\left\{\begin{array}{l} \boldsymbol{v}_{0}=\boldsymbol{u}_{0} \neq \mathbf{0} \\ \boldsymbol{v}_{k}=\boldsymbol{A} \boldsymbol{u}_{k-1},\\ \boldsymbol{u}_{k}=\frac{\boldsymbol{v}_{k}}{\max \left(\boldsymbol{v}_{k}\right)} \end{array} \quad(k=1,2,\cdots),\right.$$

有

$$\lim _{k \rightarrow \infty} \boldsymbol{u}_{k}=\frac{\boldsymbol{x}_{1}}{\max \left(\boldsymbol{x}_{1}\right)},\quad \lim _{k \rightarrow \infty} \max \left(\boldsymbol{v}_{k}\right)=\lambda_{1}$$

原点平移法加速

应用幂法计算 $\boldsymbol{A}$ 的主特征值时，其收敛速度主要由比值 $r=\frac{\lambda_{1}}{\lambda_{2}}$ 来决定，当 $r$ 接近于 $1$ 时，收敛可能很慢

引进矩阵 $\boldsymbol{B}=\boldsymbol{A}-p \boldsymbol{I}$，其中 $p$ 为选择参数。设 $\boldsymbol{A}$ 的特征值为 $\lambda_{1},\lambda_{2},\cdots,\lambda_{n}$，则 $\boldsymbol{B}$ 的相应特征值为 $\lambda_{1}-p,\lambda_{2}-p,\cdots,\lambda_{n}-p$，而且 $\boldsymbol{A}$，$\boldsymbol{B}$ 的特征向量相同。

如果需要计算 $\boldsymbol{A}$ 的主特征值 $\lambda_{1}$，就要选择适当的 $p$ 使 $\lambda_{1}-p$ 仍然是 $\boldsymbol{B}$ 的主特征值，且使

$$\left|\frac{\lambda_{2}-p}{\lambda_{1}-p}\right|<\left|\frac{\lambda_{2}}{\lambda_{1}}\right|$$

对 $\boldsymbol{B}$ 应用幂法，使得在计算 $\boldsymbol{B}$ 的主特征值 $\lambda_{1}-p$ 的过程中得到加速。这种方法通常称为原点平移法

当希望计算 $\lambda_{n}$ 时，应选择 $p=\frac{\lambda_{1}+\lambda_{n-1}}{2}=p^{*}$，使得应用幂法计算 $\lambda_{n}$ 得到加速

Rayleigh 商加速法

对称矩阵 $\boldsymbol{A}$ 的 $\lambda_{1}$ 及 $\lambda_{n}$ 可用 Rayleigh 商的极值来表示。

设 $\boldsymbol{A} \in \mathbf{R}^{n \times n}$ 为对称阵，特征值满足 $\left|\lambda_{1}\right|>\left|\lambda_{2}\right| \geqslant\left|\lambda_{3}\right| \geqslant \cdots \geqslant\left|\lambda_{n}\right|$，对应的特征向量满足 $\left(\boldsymbol{x}_{\boldsymbol{i}},\boldsymbol{x}_{j}\right)=\delta_{i j}$，应用幂法计算 $\boldsymbol{A}$ 的主特征值 $\lambda_{1}$，则规范化向量 $\boldsymbol{u}_{k}$ 的 Rayleigh 商给出 $\lambda_{1}$ 的较好的近似，即

$$\frac{\left(\boldsymbol{A} \boldsymbol{u}_{k},\boldsymbol{u}_{k}\right)}{\left(\boldsymbol{u}_{k},\boldsymbol{u}_{k}\right)}=\lambda_{1}+O\left(\left(\frac{\lambda_{2}}{\lambda_{1}}\right)^{2 k}\right)$$

反幂法

计算 $\boldsymbol{A}$ 的按模最小的特征值 $\lambda_{n}$ 的问题就是计算 $\boldsymbol{A}^{-1}$ 的按模最大的特征值问题。对 $\boldsymbol{A}^{-1}$ 应用幂法迭代法（称为反幂法），可求得矩阵 $\boldsymbol{A}^{-1}$ 的主特征值 $1 / \lambda_{n}$，从而求得 $\boldsymbol{A}$ 的按模最小的特征值 $\lambda_{n}$

反幂法迭代公式为任取初始向量 $v_{0}=\boldsymbol{u}_{0} \neq \mathbf{0}$，构造向量序列

$$\left\{\begin{array}{l} \boldsymbol{v}_{k}=\boldsymbol{A}^{-1} \boldsymbol{u}_{k-1},\\ \boldsymbol{u}_{k}=\frac{\boldsymbol{v}_{k}}{\max \left(\boldsymbol{v}_{k}\right)} \end{array} \quad(k=1,2,\cdots) \right.$$

迭代向量 $\boldsymbol{v}_{k}$ 可以通过解方程组 $\boldsymbol{A} \boldsymbol{v}_{k}=\boldsymbol{u}_{k-1}$ 求得。

对任何初始非零向量 $\boldsymbol{u}_{0}=\boldsymbol{v}_{0} \quad\left(a_{n} \neq 0\right)$，由反幂法构造的向量序列 $\left\{\boldsymbol{v}_{k}\right\},\left\{\boldsymbol{u}_{k}\right\}$ 满足

$$\lim _{k \rightarrow \infty} \boldsymbol{u}_{k}=\frac{\boldsymbol{x}_{n}}{\max \left(\boldsymbol{x}_{n}\right)},\quad \lim_{k \rightarrow \infty} \max \left(\boldsymbol{v}_{k}\right)=\frac{1}{\lambda_{n}}$$

收敛速度的比值为 $\left|\frac{\lambda_{n}}{\lambda_{n-1}}\right|$

原点平移加速

设 $\boldsymbol{A} \in \mathbf{R}^{n \times n}$ 有 $n$ 个线性无关的特征向量，$\boldsymbol{A}$ 的特征值及对应的特征向量记作 $\lambda_{i}$ 及 $\boldsymbol{x}_{i}(i=1,2,\cdots,n)$。$p$ 为 $\lambda_{j}$ 的近似值，$(\boldsymbol{A}-p \boldsymbol{I})^{-1}$ 存在，且 $\left|\lambda_{j}-p\right|<\left|\lambda_{i}-p\right|(i \neq j)$，$\boldsymbol{u}_{0}=\sum_{i=1}^{n} a_{i} x_{i} \neq 0$ 为给定的初始向量 $\left(a_{i} \neq 0\right)$，则由反幂法迭代公式构造的向量序列 $\left\{\boldsymbol{v}_{k}\right\}$，$\left\{\boldsymbol{u}_{k}\right\}$ 满足

$$\begin{array}{c} \lim _{k \rightarrow \infty} \boldsymbol{u}_{k}=\frac{\boldsymbol{x}_{j}}{\max \left(\boldsymbol{x}_{j}\right)},\\ \lim _{k \rightarrow \infty} \max \left(\boldsymbol{v}_{k}\right)=\frac{1}{\lambda_{j}-p},\quad \mathrm{即} p+\frac{1}{\max \left(\boldsymbol{v}_{k}\right)} \rightarrow \lambda_{j} \quad(\mathrm{当} k \rightarrow \infty) \end{array}$$

且收敛速度由比值 $r=\max _{i \neq j}\left|\frac{\lambda_{j}-p}{\lambda_{i}-p}\right|$ 确定

• 选取合适的 $p$，可以求 $p$ 值附近的特征值

通常不会去求 $\boldsymbol{A}^{-1}$（计算量 $n^4$），反幂法迭代公式中的 $\boldsymbol{v}_{k}$ 是通过解方程组 $(\boldsymbol{A}-p \boldsymbol{I}) \boldsymbol{v}_{k}=\boldsymbol{u}_{k-1}$ 求得的（计算量 $n^3$）。为了节省工作量，可以先将 $(\boldsymbol{A}-p \boldsymbol{I})$ 进行三角分解，即

$$\boldsymbol{P}(\boldsymbol{A}-p \boldsymbol{I})=\boldsymbol{L} \boldsymbol{U},$$

其中 $\boldsymbol{P}$ 为某个置换矩阵，于是求 $\boldsymbol{v}_{k}$ 相当于解两个三角形方程组 $\boldsymbol{L} \boldsymbol{y}_{k}=\boldsymbol{P} \boldsymbol{u}_{k-1}$，$\boldsymbol{U} \boldsymbol{v}_{k}=\boldsymbol{y}_{k}$。
反幂法迭代公式可写为

$$\left\{\begin{array}{l} \boldsymbol{L} \boldsymbol{y}_{k}=\boldsymbol{P} \boldsymbol{u}_{k-1},\\ \boldsymbol{U} \boldsymbol{v}_{k}=\boldsymbol{y}_{k},\quad(k=1,2,\cdots)\\ \boldsymbol{u}_{k}=\frac{\boldsymbol{v}_{k}}{\max \left(\boldsymbol{v}_{k}\right)} \end{array}\right.$$

按下述方法选择 $\boldsymbol{v}_{0}=\boldsymbol{u}_{0}$ 是较好的：选 $\boldsymbol{u}_{0}$ 使

$$\boldsymbol{U} \boldsymbol{v}_{1}=\boldsymbol{L}^{-1} \boldsymbol{P} \boldsymbol{u}_{0}=(1,1,\cdots,1)^{\mathrm{T}}$$

Householder 方法

设 $\boldsymbol{A} \in \mathbf{R}^{n \times n}$，则存在一个正交阵 $\boldsymbol{R}$，使

$$\boldsymbol{R}^{\mathrm{T}} \boldsymbol{A} \boldsymbol{R}=\left(\begin{array}{cccc} T_{11} & T_{12} & \cdots & T_{1 s} \\ & T_{22} & \cdots & T_{2 s} \\ & & \ddots & \vdots \\ & & & T_{ss} \end{array}\right),$$

其中对角块为一阶或二阶矩阵，每一个一阶对角块即为 $\boldsymbol{A}$ 的实特征值，每一个二阶对角块的两个特征值是 $\boldsymbol{A}$ 的一对共轭复特征值。

一方阵 $\boldsymbol{B}$，如果当 $i>j+1$ 时有 $b_{i j}=0$，则称 $\boldsymbol{B}$ 为上 Hessenberg 阵，即

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

设向量 $\boldsymbol{w}$ 满足 $\|\boldsymbol{w}\|_{2}=1$，矩阵 $\boldsymbol{H}=\boldsymbol{I}-2 \boldsymbol{w} \boldsymbol{w}^{\mathrm{T}}$ 称为初等反射阵，记作 $H(w)$，即

$$\boldsymbol{H}(\boldsymbol{w})=\left(\begin{array}{cccc} 1-2 w_{1}^{2} & -2 w_{1} w_{2} & \cdots & -2 w_{1} w_{n} \\ -2 w_{2} w_{1} & 1-2 w_{2}^{2} & \ddots & \vdots \\ \vdots & \ddots & \ddots & -2 w_{n-1} w_{n} \\ -2 w_{n} w_{1} & \cdots & -2 w_{n} w_{n-1} & 1-2 w_{n}^{2} \end{array}\right)$$

• 初等反射阵 $\boldsymbol{H}$ 是对称阵（$\boldsymbol{H}^{\mathrm{T}}=\boldsymbol{H}$），正交阵（$\boldsymbol{H}^{\mathrm{T}} \boldsymbol{H}=\boldsymbol{I}$）和对合阵（$\boldsymbol{H}^{2}=I$），有 $n-1$ 个特征值 $1$ 和一个特征值 $-1$，特征值 $-1$ 对应的特征向量就是 $\boldsymbol{w}$
• 设 $\boldsymbol{x},\boldsymbol{y}$ 为两个不相等的 $n$ 维向量，$\|\boldsymbol{x}\|_{2}=\|\boldsymbol{y}\|_{2}$，则存在一个初等反射阵 $\boldsymbol{H}$，使 $\boldsymbol{H x}=\boldsymbol{y}$
  • 易知，$w=\frac{x-y}{\|x-y\|_{2}}$ 是使 $\boldsymbol{H x}=\boldsymbol{y}$ 成立的唯一长度等于 $1$ 的向量（不计符号）
• 设向量 $\boldsymbol{x} \in \mathbf{R}^{n}(\boldsymbol{x} \neq \mathbf{0})$，$\sigma= \pm\|\boldsymbol{x}\|_{2}$，且 $\boldsymbol{x} \neq-\sigma \boldsymbol{e}_{1}$，则存在一个初等反射阵 $$\boldsymbol{H}=\boldsymbol{I}-2 \frac{\boldsymbol{u} \boldsymbol{u}^{\mathrm{T}}}{\|\boldsymbol{u}\|_{2}^{2}} \equiv \boldsymbol{I}-\rho^{-1} \boldsymbol{u} \boldsymbol{u}^{\mathrm{T}}$$ 使 $\boldsymbol{H} \boldsymbol{x}=-\boldsymbol{\sigma} \boldsymbol{e}_{1}$，其中 $\boldsymbol{u}=\boldsymbol{x}+\sigma \boldsymbol{e}_{1}$，$\rho=\|\boldsymbol{u}\|_{2}^{2} / 2$。

QR 分解

QR 分解加速方程求解

对于方程 $\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}$，若能把 $\boldsymbol{A}$ 分解成一个正交阵 $\boldsymbol{Q}$ 和一个上三角阵 $\boldsymbol{R}$ 的乘积 $\boldsymbol{A}=\boldsymbol{Q}\boldsymbol{R}$，则方程化为 $\boldsymbol{Q}\boldsymbol{R}\boldsymbol{x}=\boldsymbol{b}$，令 $\boldsymbol{y}=\boldsymbol{R}\boldsymbol{x}$，则

$$\left\{\begin{array}{c} \boldsymbol{Q}\boldsymbol{y}=\boldsymbol{b} \\ \boldsymbol{y}=\boldsymbol{R}\boldsymbol{x} \end{array}\right.$$

第一个方程中由于 $\boldsymbol{Q}^{-1}=\boldsymbol{Q}^{\mathrm{T}}$，故 $\boldsymbol{y}=\boldsymbol{Q}^{\mathrm{T}}\boldsymbol{b}$，第二个方程直接代入可以解得 $\boldsymbol{x}$

QR 分解的方法

下面考虑用初等反射阵来正交相似约化一般矩阵和对称矩阵。设

$$\boldsymbol{A}=\left(\begin{array}{c:ccc} 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)$$

选择初等反射阵 $\boldsymbol{H}_{1}$ 把矩阵的第一列约化

$$\left\{\begin{array}{l} \sigma_{1}=\operatorname{sgn}\left(a_{11}\right)\sqrt{\sum_{i=1}^{n} a_{i 1}^{2}} \\ \boldsymbol{u}_{1}=(a_{11} + \sigma_{1}, a_{21}, \cdots, a_{n1})^{\mathrm{T}} \\ \rho_{1}=\frac{1}{2}\left\|\boldsymbol{u}_{1}\right\|_{2}^{2}=\sigma_{1}\left(\sigma_{1}+a_{11}\right) \\ \boldsymbol{H}_{1}=\boldsymbol{I}-\rho_{1}^{-1} \boldsymbol{u}_{1} \boldsymbol{u}_{1}^{\mathrm{T}} \end{array}\right.$$

于是

$$\boldsymbol{A}_{2} = \boldsymbol{H}_{1}\boldsymbol{A}=\left(\begin{array}{c:ccc} -\sigma_{1} & a^{(2)}_{12} & \cdots & a^{(2)}_{1 n} \\ 0 & a^{(2)}_{22} & \cdots & a^{(2)}_{2 n} \\ \vdots & \vdots & & \vdots \\ 0 & a^{(2)}_{n 2} & \cdots & a^{(2)}_{n n} \end{array}\right)$$

下一步，取 $\boldsymbol{A}_{2}$ 右下角的子矩阵

$$\left(\begin{array}{c:cc} a^{(2)}_{22} & \cdots & a^{(2)}_{2 n} \\ \vdots & & \vdots \\ a^{(2)}_{n 2} & \cdots & a^{(2)}_{n n} \end{array}\right)$$

计算矩阵 $\hat{\boldsymbol{H}}_{2}$ 约化这个的第一列，为了把它作用在原矩阵上，只需把左上角填上单位阵即可

$$\boldsymbol{H}_{2}= \left(\begin{array}{cc} \boldsymbol{I}_{1 \times 1} & \boldsymbol{O} \\ \boldsymbol{O} & \hat{\boldsymbol{H}}_{2} \end{array}\right)$$

从合起来可以得到

$$\boldsymbol{A}_{3}=\boldsymbol{H}_{2}\boldsymbol{H}_{1}\boldsymbol{A}=\left(\begin{array}{cc:cc} -\sigma_{1} & a^{(3)}_{12} & \cdots & a^{(3)}_{1 n} \\ 0 & -\sigma_{2} & \cdots & a^{(3)}_{2 n} \\ \vdots & \vdots & & \vdots \\ 0 & 0 & \cdots & a^{(3)}_{n n} \end{array}\right)$$

依次计算，最终可以得到 $\boldsymbol{A}=\boldsymbol{Q}\boldsymbol{R}$，其中 $\boldsymbol{Q}^{\mathrm{T}} = \boldsymbol{H}_{n-1}\boldsymbol{H}_{n-2}\cdots\boldsymbol{H}_{1}$ 是一个正交阵，$\boldsymbol{R}$ 是一个上三角阵。

• QR 分解永远存在（只要 $\boldsymbol{A}$ 是方阵）
• QR 分解唯一的条件是 $\boldsymbol{A}$ 满秩，且主对角元符号要唯一

QR 算法

设 $\boldsymbol{A}=\boldsymbol{A}_{1}=\left(a_{i j}\right) \in \mathbf{R}^{n \times n}$，且对 $\boldsymbol{A}$ 进行 $Q R$ 分解，即 $\boldsymbol{A}=\boldsymbol{Q R}$，其中 $\boldsymbol{R}$ 为上三角阵，$\boldsymbol{Q}$ 为正交阵，于是可得到一新矩阵

$$\boldsymbol{B}=\boldsymbol{R} \boldsymbol{Q}=\boldsymbol{Q}^{\mathrm{T}} \boldsymbol{A} \boldsymbol{Q}$$

显然，$\boldsymbol{B}$ 是由 $\boldsymbol{A}$ 经过正交相似变换得到，因此 $\boldsymbol{B}$ 与 $\boldsymbol{A}$ 特征值相同。再对 $\boldsymbol{B}$ 进行 QR 分解，又可得一新的矩阵，重复这过程可得到矩阵序列。

$$\begin{array}{l} \boldsymbol{A}=\boldsymbol{A}_{1}=\boldsymbol{Q}_{1} \boldsymbol{R}_{1} \\ \boldsymbol{A}_{2}=\boldsymbol{R}_{1} \boldsymbol{Q}_{1}=\boldsymbol{Q}_{1}^{\mathrm{T}} \boldsymbol{A}_{1} \boldsymbol{Q}_{1} \\ \cdots \\ \boldsymbol{A}_{k}=\boldsymbol{Q}_{k} \boldsymbol{R}_{k} \\ \boldsymbol{A}_{k+1}=\boldsymbol{R}_{k} \boldsymbol{Q}_{k}=\boldsymbol{Q}_{k}^{\mathrm{T}} \boldsymbol{A}_{k} \boldsymbol{Q}_{k} \end{array}$$

按上述递推法则构造矩阵序列 $\left\{\boldsymbol{A}_{k}\right\}$ 的过程，只要 $\boldsymbol{A}$ 为非奇异矩阵，则由 QR 算法就完全确定 $\left\{\boldsymbol{A}_{k}\right\}$。

（基本 QR 方法）设 $\boldsymbol{A}=\boldsymbol{A}_{1} \in \mathbf{R}^{n \times n}$，QR 算法为

$$\left\{\begin{array}{l} \boldsymbol{A}_{k}=\boldsymbol{Q}_{k} \boldsymbol{R}_{k} \quad\left(\boldsymbol{Q}_{k}^{\mathrm{T}} \boldsymbol{Q}_{k}=\boldsymbol{I},\boldsymbol{R}_{k} \mathrm{为上三角阵}\right),\\ \boldsymbol{A}_{k+1}=\boldsymbol{R}_{k} \boldsymbol{Q}_{k} \quad(k=1,2,\cdots), \end{array}\right.$$

且记 $\widetilde{\boldsymbol{Q}}_{k} \equiv \boldsymbol{Q}_{1} \boldsymbol{Q}_{2} \cdots \boldsymbol{Q}_{k}$，$\widetilde{\boldsymbol{R}}_{k} \equiv \boldsymbol{R}_{k} \cdots \boldsymbol{R}_{2} \boldsymbol{R}_{1}$，则有

• $\boldsymbol{A}_{k+1}$ 相似于 $\boldsymbol{A}_{k}$，即 $\boldsymbol{A}_{k+1}=\boldsymbol{Q}_{k}^{\mathrm{T}} \boldsymbol{A}_{k} \boldsymbol{Q}_{k}$
• $\boldsymbol{A}_{k+1}=\left(\boldsymbol{Q}_{1} Q_{2} \cdots \boldsymbol{Q}_{k}\right)^{\mathrm{T}} \boldsymbol{A}_{1}\left(\boldsymbol{Q}_{1} Q_{2} \cdots \boldsymbol{Q}_{k}\right)=\widetilde{\boldsymbol{Q}}_{k}^{\mathrm{T}} \boldsymbol{A}_{1} \widetilde{\boldsymbol{Q}}_{k}$
• $\boldsymbol{A}^{k}$ 的 QR 分解式为 $\boldsymbol{A}^{k}=\widetilde{\boldsymbol{Q}}_{k} \widetilde{\boldsymbol{R}}_{k}$

设 $\boldsymbol{M}_{k}=\boldsymbol{Q}_{k} \boldsymbol{R}_{k}$，其中 $\boldsymbol{Q}_{k}$ 为正交阵，$\boldsymbol{R}_{k}$ 为具有正对角元素的上三角阵，如果 $\boldsymbol{M}_{k} \rightarrow \boldsymbol{I}(k \rightarrow \infty)$，则 $\boldsymbol{Q}_{k} \rightarrow \boldsymbol{I}$，及 $\boldsymbol{R}_{k} \rightarrow \boldsymbol{I}(k \rightarrow \infty)$。

收敛性

设 $\boldsymbol{A}=\left(a_{i j}\right) \in \mathbf{R}^{n \times n}$，

• 如果 $\boldsymbol{A}$ 的特征值满足：$\left|\lambda_{1}\right|>\left|\lambda_{2}\right|>\cdots>\left|\lambda_{n}\right|>0$
• $\boldsymbol{A}$ 有标准形 $\boldsymbol{A}=\boldsymbol{X} \boldsymbol{D} \boldsymbol{X}^{-1}$，其中 $\boldsymbol{D}=\operatorname{diag}\left(\lambda_{1},\lambda_{2},\cdots,\lambda_{n}\right)$，且设 $\boldsymbol{X}^{-1}$ 有三角分解 $\boldsymbol{X}^{-1}=\boldsymbol{L} \boldsymbol{U}$（$\boldsymbol{L}$ 为单位下三角阵，$\boldsymbol{U}$ 为上三角阵）

则由 QR 算法产生的 $\left\{\boldsymbol{A}_{k}\right\}$ 本质上收敛于上三角阵，即

$$\boldsymbol{A}_{k} \xrightarrow{\mathrm{本质上}} \boldsymbol{R}=\left(\begin{array}{cccc} \lambda_{1} & \times & \cdots & \times \\ & \lambda_{2} & \ddots & \vdots \\ & & \ddots & \times \\ & & & \lambda_{n} \end{array}\right) \quad(k \rightarrow \infty)$$

• 如果对称阵 $\boldsymbol{A}$ 满足收敛性条件，则由 QR 算法产生的 $\left\{\boldsymbol{A}_{k}\right\}$ 收敛于对角阵
• 设 $\boldsymbol{A} \in \mathbf{R}^{n \times n}$，如果 $\boldsymbol{A}$ 的等模特征值中只有实重特征值或多重复的共轭特征值，则由 QR算法产生 $\left\{\boldsymbol{A}_{k}\right\}$ 本质收敛于分块上三角阵（对角块为一阶和二阶子块），且对角块每一个 $2 \times 2$ 子块给出 $\boldsymbol{A}$ 的一对共轭复特征值，每一个 $1 \times 1$ 子块给出 $\boldsymbol{A}$ 的实特征值