【数值分析】学习笔记（包含精简知识点和例题，知识点略不全面，适合考前突击）

本文链接：https://blog.csdn.net/Naaameless/article/details/117294872

数值分析

这门课是建模和研究的基础，真正学习主要是了解其原理与实现方法，需要的时候再去认真学习。

于是这里首先记录考试的公式定理，顺便记录了一些算法的原理和实现，方便查阅但不是很详细，一定要结合书本看。

参考书籍：《数值分析与算法》，清华大学出版社，喻文健；
参考网络课：《数值分析》东北大学，邵新慧；

当然可以根据自己的喜好选择书籍，和网络课，大同小异。

数值分析是啥？

数值分析是研究科学计算中各种数学问题求解的数值计算方法。
说白了，就是能运用在计算机上的各种数学问题的解法、算法。

误差、精度与计算机计算

误差分类

模型误差：对实际问题建模产生的误差
观测误差：测量误差

上面2个不可避免，非重点。

截断误差（方法误差）：采用的计算方法所产生的误差
比如 通过四舍五入得到的近似值 ，采用无限项公式的某几项计算等
舍入误差：计算机只能对有限位的数进行计算
比如表示 $\pi$ 、 $\sqrt2$ 等无限位数字时，不能完全相等

这2个是本科目研究的重点。

误差的概念

一般，设 $x$ 为精确值 $x^*$ 的一个近似值，记：
$e = x^*-x$
e 为近似值 $x$ 的绝对误差，简称误差。
如果， $\leqslant \epsilon$ ，称 $\epsilon$ 为近似值 $x$ 的绝对误差限。

就是说，
$x-\epsilon \leqslant x^* \leqslant x + \epsilon$
但是，它并不能真正反正出误差的影响，于是令： $e_r = \frac \epsilon {x^*} = \frac {x^*-x} {x^*}$
称 $e_r$ 为近似值 $x$ 的相对误差，由于 $X^*$ 一般未知，所以取值为：
$e_r = \frac \epsilon {x} = \frac {x^*-x} {x}$
于是它也有一个误差限， $\epsilon_r = \frac \epsilon {|x|}$ 称为近似值 $x$ 的相对误差限。 $|e_r| \leqslant \epsilon_r$ 。

**例：**设 $x = 1.24$ 是有精确值 $x^*$ 经过四舍五入得到的近似值，
求 $x$ 的绝对误差限和相对误差限。
解：已知， $\leqslant x^* \leqslant 1.245$
所以， $\epsilon = \frac {1.245-1.235} 2 = 0.005$ ， $\epsilon_r = \frac {0.005} {1.245} \approx 0.4\%$

定理：一般地，凡是由精确值经过四舍五入得到的近似值，其绝对误差限等于该值末位的半个单位。
如上题 $x$ 的末位为百分位，即0.01位，所以其 $\epsilon = 0.005$ 是它的一半。

有效数字与误差的关系

定义：设数 x 是数 x 的近似值，如果 x 的绝对误差限是他的某一数位的半个单位*，并且从 x 左起第一个非零数字到该数字共有 n 位，则称这 n 个数字为 x 的有效数字，也称 x 近似 X* 时具有 n 位有效数字。

要知道，任何数可写为标准浮点数形式：
$\pm 0.a_1a_2 \dots a_k \times 10^m$
m 是整数， $a_k$ 是 0 - 9 的数字， $a_1 \ne 0$ 。

那么，有效数字与绝对误差与绝对误差限的关系：
$|x^*-x| = e \leqslant \epsilon \leqslant \frac 12 \times 10^{m-n}$
m 是标准浮点数形式中的 m，n 是有效数字个数。

例：为使 $x^* = \sqrt2$ 的近似值的绝对误差小于 $10^{-5}$ ，
问应取几位有效数字？
解：近似值 x 可写为：
$\pm 0.a_1a_2 \dots a_k \times 10$ ， $a_1 \ne 0$
令 $\frac 12 \times 10^{1-n} \leqslant 10^{-5}$
故取 $n = 6$ ，即要取6位有效数字
思路：已知绝对误差限 $\epsilon < 10^{-5}$ ，由 $\epsilon = |x^*-x| \leqslant \frac 12 \times 10^{m-n}$ ，可知 $\frac 12 \times 10^{1-n} \leqslant 10^{-5}$ ，于是便可求得 n 。

例： $\pi$ 的近似值 3.14，3.141，3.15，
求他们的精确程度；他们每个位上的数字都起作用吗？
解：各近似值绝对误差为： $|\pi - 3.14 |= 0.0015...， |\pi - 3.141| = 0.0005... ， |pi - 3.15| = 0.0085...$
所以，3.141最精确，3.15误差最大
各近似值分别写为： $3.14 \times 10 ， 3.141 = 0.3141 \times 10 ，3.15 = 0.315 \times 10$
由 $|x^*-x| \leqslant \frac 12 \times 10^{m-n}$ ，
对3.14来说， $\leqslant \frac1 2 \times 10^{1-n}$ 应取 n = 3
对3.141来说， $\leqslant \frac1 2 \times 10^{1-n}$ 应取 n = 3
对3.15来说， $\leqslant \frac1 2 \times 10^{1-n}$ 应取 n = 2
故3.14、3.141有3位有效数字，3.15有2位有效数字
思路：精确程度指绝对误差，数字作用指有效数字；已知绝对误差，利用定义和公式，可求得有效位数 n 与对应的绝对误差限。

**拓展：**一般 m 已知，由公式可以看出，再知道 e $\epsilon$ n 中的一个，就能求出其它。

有效数字与相对误差与相对误差限的关系：
若 x 有 n 位有效数字，则其相对误差限为： $\frac {|x^*-x|} x = e_r \leqslant \epsilon_r \leqslant \frac 1 {2a_1} \times 10^{1-n}$
若 x 的相对误差限为： $\epsilon_r \leqslant \frac 1 {2(a_1+1)} \times 10^{1-n}$ ，则至少有 n 位有效数字。

敏感性（病态性）与数值稳定性

问题的敏感性（病态性）：指输入数据的扰动对问题解的影响程度的大小。它是问题的属性，与算法无关。
如果解的相对变化远远超过输入数据的变化，则称这个问题是 敏感的 或者 病态的 ；反之，称这个问题是 不敏感的 或者 良态的 。

量化问题敏感程度的条件数，被定义为 $\rm {cond} = \frac {||问题的解的相对变化量||} {||输入数据的相对变化量||} $ ，其中 $||\cdot||$ 表示范数

算法的稳定性：指在计算过程中的扰动对问题解的影响程度。它是算法的属性，与问题无关。
在计算过程中，舍入误差或者小扰动不被放大或者放大不严重，则称算法是稳定的（具有数值稳定性）；否则，称不稳定的（不具有算法稳定性）。

数值计算中的原则

设计算法的时候，为了避免放大误差与计算稳定，遵循以下原则：

避免两个相近的数相减：会造成相对误差变大，有效数字减小
- 解决策略：转换计算方法，如取对数、数学变换、利用其他公式，或计算机使用双倍字长运算
在求和或差时由小到大运算：防止在数量级差距过大时，因计算机精度不够的计算方式舍去了小的数值
避免绝对值太小的数作除数：可能会导致商的绝对误差放大
注意化简计算程序：减少计算量
选择数值稳定性好的算法：算法在计算舍入误差积累是可控制的，则称为数值稳定的

线性方程组的直接解法

线性方程组可以写为：

$\mathbf A\mathbf x = \mathbf b$
A 是系数矩阵，x 是 n 维向量（解向量，变量），b 是 n 维向量（已知向量）。

设 $\mathbf A = a_{ij}$ ， $i=1,\dots, m$ , $1,\dots,n$ ，
那么记在 $a_{ij}$ 位置上，第 k 次出现的数值为 $a_{ij}^{(k)}$ ，同理有 $b_{i}^{(k)}$ 。

当 m > n，称 超定方程组 ，一般无解，但可求出最小二乘解；
当 m < n，方程组一般有无穷多解，实际应用中常和约束条件一起构成约束优化问题；
当 m = n，方程组一般有唯一解，是本课程主要研究的情况。

直接解法适用于：A 为低阶稠密矩阵的情况（低阶：阶数 n 较小；稠密：零元素较少）
其中较为有效的是 列主元消去法 。

顺序高斯(Gauss)消去法

原理：对增广矩阵 $(\mathbf A,\mathbf b)$ 进行初等行变换，直到化成 阶梯矩阵 ，然后回代计算出解。

思路：假设在消去过程中主对角元素始终不为 0 ；
不断利用每行主对角元素消去下面同列的元素，逐渐使矩阵变成阶梯矩阵（上三角矩阵）；
完成后回代求解。

优点：计算量较小，结果准确

缺点：主对角元素在消去过程中不能出现 0 的情况 $\Rightarrow$ 矩阵A的各界顺序主子式都不为零；
可能出现小数（主对角元素）作为除数，或者大数与小数加减的情况，即可能数值稳定性弱

列主元高斯消去法

原理：利用选择主元的技术，改善顺序高斯消去法。

思路：假设 A 为非奇异矩阵，在进行顺序高斯消元之前，先在当前主元所在列，从当前行及其下方的元素中选出最大的元素（ $|a_{kk}^{(k)}| = \max \limits_{k \le i \le n} |a_{ik}^{(k)}|$ ），交换当前行与最大元素所在行的位置，从而确保主元不等于零，且不会太小。

优点：让高斯消元法遇到主元素为 0 的时候能够继续；
同时也为了减轻当主元是小数时，作为除数带来的误差放大;
前提条件比顺序高斯消去法宽松。

缺点：要求 $|\mathbf A| \ne 0$ 。

矩阵三角分解法（LU分解）

原理（LU分步分解法）

在高斯消元法中，一直使用 初等行变换 进行消元，而根据线性代数的知识我知道，初等行变换等价于左乘一系列 初等矩阵 。于是可利用其将原矩阵 A 分解：
$\begin{aligned} &\mathbf A^{(n)} = \mathbf L_{n-1}\mathbf L_{n-2} \dots \mathbf L_1A^{(1)};\\ \Rightarrow& \mathbf A^{(1)} = \mathbf L_1^{-1}\mathbf L_2^{-1} \dots \mathbf L_{n-1}^{-1} \mathbf A^{(n)}\\ \Rightarrow& \mathbf A = \mathbf L\mathbf U \end{aligned}$
其中，左乘的 单位下三角矩阵（主对角线为1，上三角为0；单位上三角矩阵同理），都能被求出来。其中第 k 个初等矩阵为：
$\mathbf L_k = \begin{bmatrix} 1 & \quad & \quad & \quad & \quad & \quad \\ \quad & \ddots & \quad & \quad & \quad & \quad \\ \quad & \quad & 1 & \quad & \quad & \quad \\ \quad & \quad & -l_{k+1k} & 1 & \quad & \quad \\ \quad & \quad & \vdots & \quad & \ddots & \quad \\ \quad & \quad & -l_{nk} & \quad & \quad & 1 \\ \end{bmatrix} \quad ,k = 1,2,\dots,n-1$
于是，
$\mathbf L =\mathbf L_1^{-1}\mathbf L_2^{-1} \dots \mathbf L_{n-1}^{-1} \\ \mathbf U = \mathbf A^{(n)}$
可见，矩阵 L 是个 单位下三角矩阵 ，U 是个 上三角矩阵 。

直接LU分解法（三角分解，Doolittle分解法）

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

故，可以放心，对于各阶顺序主子式不为零的n阶方阵A，分步求L、U矩阵与 直接求出L、U矩阵结果相同，后者更快捷。

思路：

计算顺序：先计算 U 矩阵的第一行，在计算 L 矩阵的第一列，再计算 U 的一行，再计算 L 的一列…直到计算到 $u_{nn}$
在计算LU矩阵某位置的元素时，都是根据同一位置的系数矩阵的元素与分解矩阵的乘积规则 $a_{ij} = l_iu_j$ 得到的。
数据存放的时候，可以发现 L 与 U 矩阵正好可以存在原 A 矩阵中，不需要另外开辟存储空间
回代：先回代 L 矩阵计算出 $\mathbf B^{(n)}$ ，再结合 $\mathbf B^{(n)}$ 与 U 矩阵计算出 X

对于系数矩阵 A，使用矩阵三角分解的 Doolittle 分解法，设：
$\begin{bmatrix} a_{11} & a_{12} & \dots & a_{1n} \\ a_{21} & a_{22} & \dots & a_{2n} \\ \vdots & \vdots & & \vdots \\ a_{n1} & a_{n2} & \dots & a_{nn} \\ \end{bmatrix} = \begin{bmatrix} 1 & \quad & \quad & \quad & \quad \\ l_{21} & 1 & \quad & \quad & \quad \\ l_{31} & l_{32} & 1 & \quad & \quad \\ \vdots & \vdots & \vdots & \ddots & \quad \\ l_{n1} & l_{n2} & l_{n3} & \dots & 1 \\ \end{bmatrix} \begin{bmatrix} u_{11} & u_{12} & \dots & u_{1n} \\ & u_{22} & \dots & u_{2n} \\ & & \ddots & \vdots \\ & & & u_{nn} \\ \end{bmatrix}$

则，计算方法为：

$\begin{aligned} & u_{1j} = a_{1j} \quad j = 1,2,\dots,n \\ & l_{i1} = a_{i1} \div u_{11} \quad i = 2,3,\dots,n \\ & 对 k=2,3,\dots,n，有 \\ & u_{kj} = a_{kj} - \sum \limits_{m=1}^{k-1} l_{km} u_{mj} \quad j = k,k+1,\dots,n \\ & l_{ik} = (a_{ik} - \sum \limits_{m=1}^{k-1} l_{im} u_{mk}) \div u_{kk} \quad i=k+1,k+2,\dots,n \end{aligned}$
回代：由
$\begin{bmatrix} 1 & \quad & \quad & \quad & \quad \\ l_{21} & 1 & \quad & \quad & \quad \\ l_{31} & l_{32} & 1 & \quad & \quad \\ \vdots & \vdots & \vdots & \ddots & \quad \\ l_{n1} & l_{n2} & l_{n3} & \dots & 1 \\ \end{bmatrix} \begin{bmatrix} y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_n \end{bmatrix}= \begin{bmatrix} b_1 \\ b_2 \\ b_3 \\ \vdots \\ b_n \end{bmatrix} \\ \begin{bmatrix} u_{11} & u_{12} & \dots & u_{1n} \\ & u_{22} & \dots & u_{2n} \\ & & \ddots & \vdots \\ & & & u_{nn} \\ \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix}= \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}$
得
$\begin{aligned} & y_1 = b_1 \\ & y_k = b_k - \sum \limits_{i=1}^{k-1} l_{ki} y_i \quad k = 2,3,\dots,n \\ & x_n = y_n \div u_{nn} \\ & x_i = (y_i - \sum \limits_{j=i+1}^n u_{ij} x_j) \div u_{ii} \quad i=n-1,n-2,\dots,1 \end{aligned}$

优点：它与高斯消去法的计算量基本相同，但是在等式右边的项不同的方程组求解中，可以大大节省计算量。

例：使用三角分解法求解方程组 $\left \lbrace \begin{matrix} x_1+2x_2-3x_3=1 \\ 2x_1-x_2+3x_3=5 \\ 3x_1-2x_2+2x_3=1 \end{matrix} \right.$ 。

思路：系数矩阵为 $\mathbf A = \begin{bmatrix} 1 & 2 & -3 \\ 2 & -1 & 3 \\ 3 & -2 & 2 \\ \end{bmatrix}$

首先根据 $u_{1j} = a_{1j}$ 写出 U 的第一行 $\begin{bmatrix} 1 & 2 & -3 \\ & & \\ & & \\ \end{bmatrix}$ ，再根据 $l_{i1} = a_{i1}$ 写出 L 的第一列 $\begin{bmatrix} 1 & 2 & -3 \\ 2 & & \\ 3 & & \\ \end{bmatrix}$

然后根据 $\mathbf A = \mathbf L \mathbf U$ ，先计算 U 的行
$a_{22} = l_{21}u_{12} + l_{22}u_{22} + l_{23}u_{32} \\ -1 = 2 \times 2 + 1 \times u_{22} + 0 \times 0 \\ u_{22} = -5$
同理，得 $u_{23} = 9$ 。
再计算 L 的列
$a_{32} = l_{31}u_{12} + l_{32}u_{22} + l_{33}u_{32} \\ -2 = 3 \times 2 + l_{32} \times -5 + 1 \times 0 \\ l_{32} = \frac 85$
这样 L 矩阵就求解完了， $\begin{bmatrix} 1 & 2 & -3 \\ 2 & -5 & 9 \\ 3 & \frac 85 & \\ \end{bmatrix}$

再继续计算 U 的行，得 $u_{33} = \frac {-17}5$

于是 $\begin{bmatrix} 1 & 2 & -3 \\ 2 & -5 & 9 \\ 3 & \frac 85 & \frac {-17}5 \\ \end{bmatrix}$ ，即 $\mathbf L = \begin{bmatrix} 1 & & \\ 2 & 1 & \\ 3 & \frac 85 & 1 \\ \end{bmatrix} \quad \mathbf U = \begin{bmatrix} 1 & 2 & -3 \\ & -5 & 9 \\ & & - \frac {17}5 \\ \end{bmatrix}$

开始回代，先解 $\begin{bmatrix} 1 & & \\ 2 & 1 & \\ 3 & \frac 85 & 1 \\ \end{bmatrix} \begin{bmatrix} y_1 \\ y_2 \\ y_3 \end{bmatrix} = \begin{bmatrix} 1 \\ 5 \\ 1 \end{bmatrix}$ ，得 $\begin{bmatrix} y_1 \\ y_2 \\ y_3 \end{bmatrix} = \begin{bmatrix} 1 \\ 3 \\ - \frac {34}5 \end{bmatrix}$

再解 $\begin{bmatrix} 1 & 2 & -3 \\ & -5 & 9 \\ & & -\frac {17}5 \\ \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} 1 \\ 3 \\ - \frac {34}5 \end{bmatrix}$ ，得 $\begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} 1 \\ 3 \\ 2 \end{bmatrix}$

平方根法（正定矩阵的 Cholesky 分解法）

定理：设 $\mathbf A$ 为实对称正定矩阵，则存在唯一的非奇异下三角矩阵 $\mathbf L$ ，使得 $\mathbf A = \mathbf L \mathbf {L^T}$ ，且 $\mathbf L$ 的主对角元素均大于零。

于是可将 $\mathbf A \mathbf x = \mathbf b$ 转化为 $\mathbf L \mathbf y = \mathbf b$ ， $\mathbf {L^T} \mathbf x = \mathbf y$ ，称平方根法，或 Cholesky 分解法。且 L 矩阵的主对角元素取正值。

若记 $\mathbf L = (l_{ij})$ ，则对 $k = 1, 2, . . ., n$ 有：
$\begin{aligned} & l_{kk} = (a_{kk}-\sum \limits_{m=1}^{k-1} l_{km}^2)^{\frac 12} = \sqrt {a_{kk}-\sum \limits_{m=1}^{k-1} l_{km}^2} \\ & l_{ik} = (a_{ik}-\sum \limits_{m=1}^{k-1} l_{im} l_{km}) / l_{kk} \quad i = k+1,\dots,n \end{aligned}$