研一课程笔记-数值分析

绪论与预备知识

概论

衡量数值算法的几个标准

衡量算法的优劣有两个标准，一是要有可靠的理论基础，包括收敛性、数值稳定性和算法精度等；二是要有较好的计算复杂度。

• 数值稳定性

  • 数值稳定性即算法对舍入误差的敏感性。舍入误差对计算结果的精确性影响小的算法，具有较好的数值稳定性；反之，算法的数值稳定性差。

• 算法精度

  • 数值算法的精度衡量了近似解对精确解的近似程度。算法的精度越高并不一定导致近似解的误差越小。精度的具体定义依赖于所考虑的问题类型。在偏微分方程的数值算法中我们一般考虑收敛精度，它指的是算法的收敛速率。

• 计算复杂度

  • 计算复杂度指的是通过数值计算解决问题的困难程度。最常见的是时间复杂度和空间复杂度

• 收敛性

  • 算法的收敛性一般指数值解逼近精确解的性质，而对于不同的问题，其定义也是有区别的。例如在迭代算法中，收敛性指的是近似解能够随着迭代过程趋向于精确解，它通常与迭代矩阵的谱半径相关。在偏微分方程数值算法中，收敛性指的是随网格尺寸等参数变化的数值解的收敛性质，它与算法的收敛精度密切相关。

误差理论

误差按照它们的来源可分为以下四种

1. 模型误差：反映实际问题有关量之间关系的计算公式，即数学模型，通常只是近似的。由此产生的数学模型的解与实际问题的解之间的误差称为模型误差。

2. 观测误差：数学模型中包含的某些参数（如时间、长度、电压等等）往往通过观测而获得。由观测得到的数据与实际的数据之间是有误差的。这种误差称为观测误差。

3. 截断误差求解数学模型所用的数值计算方法如果是一种近似的方法，那么只能得到数学模型的近似解，由此产生的误差称为截继误差或方法误差。例如，由Taylor (泰勒)公式，函数 $f(x)$ 可表示为

$$f(x)=f(0)+f'(0)x+\frac{f''(0)}{2!}x^2+\cdots+\frac{f^{(n)}(0)}{n!}x^n+\frac{f^{(n+1)}(\theta x)}{(n+1)!}x^{n+1},(0<\theta<1)$$

近似即：（此近似公式的误差就是截断误差）

$$f(x)\approx f(0)+f'(0)x+\frac{f''(0)}{2!}x^2+\cdots+\frac{f^{(n)}(0)}{n!}x^n$$

4. 舍入误差：由于计算机的字长有限，参加运算的数据以及运算结果在计算机上存放会产生误差

数值计算的误差度量

对于准确值$x$、近似值 $a$ 和误差限 $\epsilon$ ：

• $e=x-a$ ，称 $e$ 为近似值 $a$ 的绝对误差，简称误差。

• 如果 $|e|$ 的一个上界已知，记为 $\epsilon$ ，即 $|e|\leq \epsilon$ 则称 $\epsilon$ 为近似值 $a$ 的绝对误差限或绝对误差界，简称误差限或误差界。

• $e_r=\frac ex=\frac{x-a}x$ ，称 $e_r$ 为近似值 $a$ 的相对误差。由于 $x$ 未知，实际上总把 $\frac e a$作为 $a$ 的相对误差，并且也记为$e_r=\frac ea=\frac{x-a}a$。一般用百分比表示

• $|e_r|$ 的上界，即 $\epsilon_r=\frac{\epsilon}{|a|}$ 称为近似值 $a$ 的相对误差限或相对误差界。

• 有效数字：由准确值经过四舍五入得到的近似值（经过对于某一绝对误差限舍入后），从它的末位数字到第一位非零数字都是有效数字

  • 严格来说：如果 $a$ 的绝对误差限是它的某一位的半个单位，并且从该位到它的第一位非零数字共有 $n$ 位，则称用 $a$ 近似 $x$ 时具有 $n$ 位有效数字。

PS. $a-\epsilon\leq x\leq a+\epsilon$，即 $x=a\pm\epsilon.$

函数值的误差估计

给定多元函数$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^*,\cdots,x_n^*\right)$

绝对误差

在点 $x=\left(x_1,x_2,\cdots,x_n\right)$ 进行泰勒展开，由

$$f(x_{1}, x_{2}, \ldots, x_{n}) = f(a_{1}, a_{2}, \ldots, a_{n}) + \sum_{i=1}^{n} \frac{\partial f}{\partial x_{i}} (x_{i} - a_{i}) + \frac{1}{2!} \sum_{i,j=1}^{n} \frac{\partial^{2} f}{\partial x_{i} \partial x_{j}} (x_{i} - a_{i})(x_{j} - a_{j}) + \cdots$$

其中， $\frac{\partial f}{\partial x_{i}}$ 表示函数 $f$ 对变量 $x_{i}$ 的偏导数，$\frac{\partial^{2} f}{\partial x_{i} \partial x_{j}}$ 表示函数$f$ 对变量 $x_{i}$ 和 $x_{j}$ 的混合偏导数，依此类推。所以：

$$\begin{aligned}A^{*}-A&=f\left(x_{1}^{*},\cdots,x_{n}^{*}\right)-f\left(x_{1},x_{2},\cdots,x_{n}\right)\\&\approx\sum_{j=1}^n\frac{\partial f(x)}{\partial x_j}\left(x_j^*-x_j\right)\end{aligned}$$

PS.这个等式也可f以用中值定理直接得出

由 $e(x_i)=x_j^*-x_j$ ，也即：

$$e(A)\approx\sum_{j=1}^n\frac{\partial f(x)}{\partial x_j}e\left(x_j\right),|e(A)|\leqslant\sum_{j=1}^n\left|\frac{\partial f(x)}{\partial x_j}\right||e\left(x_j\right)|$$

PS.

相对误差

$$e_r(A)=\frac{e(A)}{A}\approx\sum_{j=1}^n\frac{x_j}{f(x)}\frac{\partial f(x)}{\partial x_j}e_r\left(x_j\right)$$

• $e_r (x_j) = \frac{e(x_j)}{x_j}$ ，即 $e(x_j)= e_r (x_j)\cdot {x_j}$
• 因子 $\frac{x_j}{f(x)}\frac{\partial f(x)}{\partial x_j}$ 反映 $x_j$ 相对误差 $e_r (x_j)$ 对相对误差 $e_r(A)$ 影响的程度

函数值误差的代数运算

• $f(x,y)=xy$
  • $e(xy)\approx ye(x)+xe(y)$
  • $e_r(xy)\approx e_r(x)+e_r(y)$
• $f(x,y)=x/y$
  • $e(x/y)\approx\frac{1}{y}e(x)-\frac{x}{y^{2}}e(y)$
  • $e_{r}(x/y)\approx e_{r}(x)-e_{r}(y)$
• $f(x,y)=x\pm y$
  • $e(x\pm y)\approx e(x)\pm e(y)$
  • $e_{r}(x\pm y)\approx\frac{x}{x\pm y}e_{r}(x)\pm\frac{y}{x\pm y}e_{r}(y)$ ， $x\pm y\neq0$
• $f(x)=\sqrt{x}$
  • $e_r(\sqrt{x})\approx\frac{1}{2}e_r(x)$
• $f(x)=x^{n}$
  • $e_{r} (x^{n})\approx ne_{r}(x)$

几点注意：

• 对于 $e_{r}(x- y)\approx\frac{x}{x- y}e_{r}(x)-\frac{y}{x- y}e_{r}(y)$ ，如果 $x$ 与 $y$ 非常接近，会出现很大计算误差
• 同样的，对于 $e(x/y)\approx\frac{1}{y}e(x)-\frac{x}{y^{2}}e(y)$ ，若 $y\to 0$，同理

算法在数值计算中应注意的几个问题

• 要有数值稳定性。即能控制舍入误差的传播
• 两数相加要防止较小的数加不到较大的数中所引起的严重后果
• 避免两个相近的近似值相减，以免严重损失有效数字
• 除法运算中，要尽量避免除数的绝对值远小于被除数的绝对值

范数

向量范数

定义向量大小的量，又称为向量的模。是定义在 $R^n$ 上的实值函数 $\|\cdot\|$，满足：

• $\|x\|\geq 0$ 且仅 $x=0,\|x\|=0$
• $\|kx\|=|k|\cdot\|x\|,k\in R$
• $\|x+y\|\leq\|x\|+\|y\|$

证明是向量范数，只需要满足上面三个条件

三种向量范数($p-$范数)

对 $R^n$ 中的任一向量 $x=(x_{1},x_{2},\cdots,x_{n})^{\mathrm{T}}$

$$\|x\|_{1}=\sum_{i=1}^{n}|x_{i}|\\\|x\|_{2}=\sqrt{\sum_{i=1}^{n}x_{i}^{2}}\\\|x\|_{\infty}=\max_{1\leqslant i\leqslant n}|x_{j}|$$

范数等价性定理

设 $\|\cdot\|_\alpha,\|\cdot\|_\beta$ 是 $R^n$ 上的任意两种向量范数，则存在与向量 $x$ 无关的常数 $m$ 和 $M(0 < m < M)$， 使下列关系成立

$$m\|x\|_\alpha\leqslant\|x\|_\beta\leqslant M\|x\|_\alpha,\quad\forall x\in\mathbb{R}^n$$

此外：

摘自 张绍飞，赵迪. 矩阵论教程[M]. 机械工业出版社，2012. p98

矩阵范数

用于定义矩阵“大小”的量，又称为矩阵的模，定义在 $R^{n\times n}$ 上的实值函数 $\|\cdot\|$，满足对任意矩阵 $A,B$：

• $\|A\|\geq 0$ 且仅 $A=O,\|A\|=0$
• $\|kA\|=|k|\cdot\|A\|,k\in \mathbb R$
• $\|A+B\|\leq\|A\|+\|B\|$
• $\|AB\|\leq\|A\|\cdot\|B\|$

矩阵范数相容性

对于给定的向量范数 $\|\cdot\|$ 和矩阵范数 $\|\cdot\|$ ， 如果对任一个 $x\in \mathbb R^n$ 和任一个 $A\in \mathbb R^{n\times n}$， 满足 $\|Ax\|\leqslant\|A\|\|x\|$ 则称所给的矩阵范数与向量范数是相容的

当在同一个问题中需要同时使用矩阵范数和向量范数时，这两种范数应当是相容的

矩阵的算子范数定理

设在 $R^n$ 中给定了一种向量范数，对任一矩阵 $A\in\mathbb{R}^{n\times n}$， 令：

$$\|A\|=\max_{\|x\|=1}\|Ax\| \tag{1.3.1}$$

则由式1.3.1定义的 $\|\cdot\|$是一种矩阵范数，并且它与所给定的向量范数相容

称式1.3.1所定义的矩阵范数为从属于所给定向量范数的矩阵范数，又称为矩阵的算子范数

$\max_{\|x\|=1}\|Ax\|$ 表示矩阵 $A$ 作用在单位范数向量 $x$ 上所得到的结果的最大范数。具体来说：

• 对于给定的矩阵 $A$，我们考虑所有单位范数向量 $x$，即满足 $\|x\|=1$ 的向量。
• 我们将矩阵 $A$ 作用在每一个单位范数向量 $x$ 上，得到向量 $Ax$。
• 然后计算每个 $Ax$ 的范数 $\|Ax\|$。
• 最终，我们找到使得 $\|Ax\|$ 最大的单位范数向量 $x$，即求解 $\max_{\|x\|=1}\|Ax\|$。

证明为什么是一种矩阵范数：

设给定的向量范数为 $\|\cdot\|_p$，则从属于向量范数 $\|\cdot\|_p$ 的矩阵范数仍记为 $\|\cdot\|_p$，即

$$\|A\|_p=\max_{\|x\|_p=1}\|Ax\|_p$$

矩阵 $p-$范数

设 $A=[a_{ij}]\in\mathbb{R}^{n\times n}$，则

$$\begin{gathered} \|A\|_{1}=\max_{1\leqslant j\leqslant n}\sum_{i=1}^{n}|a_{ij}| \\ \|A\|_{2}=\sqrt{\lambda_{\max}\left(A^{\mathrm{T}}A\right)} \\ \|A\|_{\infty}=\max_{1\leqslant i\leqslant n}\sum_{j=1}^{n}|a_{ij}| \end{gathered}$$

其中 $\lambda_{\mathrm{max}} (A^{\mathrm{T}}A)$表示矩阵 $A^{\mathrm{T}}A$ 的最大特征值。

矩阵范数 $\|A\|_{1},\|A\|_{2},\|A\|_{\infty}$ 又分别称为矩阵的列范数、谱范数和行范数

佛罗贝尼乌斯范数

另外还有一种常用的矩阵范数，就是Frobenius(佛罗贝尼乌斯) 范数，又称为Euclid 范数，

$$\|A\|_\mathrm{F}=\sqrt{\sum_{i,j=1}^na_{ij}^2}$$

其中 $A=[a_{ij}]\in\mathbb{R}^{n\times n}$，又可以记作 $\|A\|_\mathrm{E}$

可以证明,$\|\cdot\|_\mathrm{F}$ 与向量范数 $\|\cdot\|_{2}$ 相容，即

$$\|Ax\|_2\leqslant\|A\|_\text{F}\|x\|_2,\quad A\in\mathbb{R}^{n\times n},x\in\mathbb{R}^n$$

但 $\|\cdot \|_\mathrm{F}$ 不从属于任何向量范数，即不是算子范数

算子范数的估计

单位矩阵 $I$ 的任何一种算子范数都有

$$\|I\|=\max\limits_{\|x\|=1}\|Ix\|=1$$

---

设矩阵 $A \in R^{n\times n}$ 的某种范数 $\left\|A\right\|<1$ ，则 $I\pm A$ 为非奇异矩阵，并且当该种范数为算子范数时，下式成立

$$\left\|(I\pm A)^{-1}\right\|\leqslant\frac1{1-\|A\|}$$

若 $I\pm A$ 奇异，即 $(I\pm A)\cdot x=0$，展开并移项，即 $x=\mp Ax$，取范数，即 $\|x\|=\|Ax\| \leq \|A\|\cdot \|x\|$

即 $\|A\|\geq1$，矛盾

当然这里还涉及到一个点，即矩阵范数相容性，即满足 $\|Ax\|\leqslant\|A\|\|x\|$

此外也可以通过这种方式来证明：

取

$$e_i=\left(\begin{array}{c}0\\1\\0\end{array}\right)$$

则：

$$x\cdot e_i=\begin{bmatrix} 0 & x_1 & 0 \\ 0 & x_2 & 0 \\ 0 & x_3 & 0 \end{bmatrix}$$

取范数，即 $\|x\|=\|x\cdot e_i\|$，相当于我们构建了一个向量范数与矩阵范数的等式

那么 $\|Ax\|=\|Ax\cdot e_i\|\leq\|A\|\cdot\|x\cdot e_i\|=\|A\|\cdot\|x\|$

由 $(I-A)(I-A)^{-1}=I$ ，则 $(I-A)^{-1}-A\cdot (I-A)^{-1}=I$，即 $(I-A)^{-1}=A\cdot (I-A)^{-1}+I$

两边同取算子范数：$\|(I-A)^{-1}\|=\|A\cdot (I-A)^{-1}\|+\|I\|\leq\|A\|\cdot\|(I-A)^{-1}\|+\|I\|$

移项，即 $(1-\|A\|)\|(I-A)^{-1}\|\leq\|I\|=1$，而 $\|A\|<1$，故上式得证

PS.一个方阵是奇异的，如果它的行列式为零。非奇异矩阵则是行列式非零的方阵。

1. 逆矩阵：奇异矩阵是不可逆的，因为逆矩阵的存在要求矩阵的行列式不为零。因此，奇异矩阵没有逆矩阵。

2. 秩：奇异矩阵的秩小于其阶数。一个 $n \times n$ 的奇异矩阵的秩最多为 $n-1$。

3. 零空间：奇异矩阵的零空间（也称为核）不仅包含零向量，还包含其他非零向量。这是因为奇异矩阵将某些向量映射到零向量。

4. 特征值：奇异矩阵的特征值中至少有一个为零。这是因为特征值是行列式与矩阵迹的根，而奇异矩阵的行列式为零。

5. 解的存在性：对于线性方程组 $Ax = b$，其中 $A$ 是奇异矩阵，解可能不存在，或者有无穷多个解。这取决于向量 $b$ 是否位于矩阵 $A$ 的列空间中。

线性方程组的解法

主要讨论如下形式的线性方程组的求解问题.

$$\left.\left\{\begin{array}{c}a_{11}x_1+a_{12}x_2+\cdots+a_{1n}x_n=b_1,\\a_{21}x_1+a_{22}x_2+\cdots+a_{2n}x_n=b_2,\\\vdots\\a_{n1}x_1+a_{n2}x_2+\cdots+a_{nn}x_n=b_n,\end{array}\right.\right.\tag{2.1.1}$$

也即：

$$Ax=b$$

其中

$$\left.A=\left(\begin{array}{ccc}a_{11}&\cdots&a_{1n}\\\vdots&&\vdots\\a_{n1}&\cdots&a_{nn}\end{array}\right.\right)\in\mathbb{R}^{n\times n},\quad x=\left(\begin{array}{c}x_1\\\vdots\\x_n\end{array}\right)\in\mathbb{R}^n,\quad b=\left(\begin{array}{c}b_1\\\vdots\\b_n\end{array}\right)\in\mathbb{R}^n.$$

设系数矩阵 $A$ 非奇异，即 $det(A) \neq 0$，则方程组有唯一解向量 $x$.

Gauss 消去法

由消元和回代两部分组成

消元即是对增广矩阵 $[A,b]$ 做初等行变换，经过 $n-1$ 次消元得到

$$\left.\left(A^{(n)},b^{(n)}\right)=\left(\begin{array}{ccccc}a_{11}^{(1)}&\cdots&\cdots&a_{1n}^{(1)}&b_{1}^{(1)}\\&a_{22}^{(2)}&\cdots&a_{2n}^{(2)}&b_{2}^{(2)}\\&&\ddots&\vdots&\vdots\\O&&&a_{nn}^{(n)}&b_{n}^{(n)}\end{array}\right.\right)$$

即线性方程组：

$$\begin{aligned}\begin{cases}a_{11}^{(1)}x_1+a_{12}^{(1)}x_2+\cdots+a_{1n}^{(1)}x_n=b_1^{(1)},\\a_{22}^{(2)}x_2+\cdots+a_{2n}^{(2)}x_n=b_2^{(2)},\\\vdots\\a_{nn}^{(n)}x_n=b_n^{(n)}.\end{cases}\end{aligned}$$

通过逐步回代，依次求出 $x_{n-1},x_{n-2},\ldots,x_1$

顺序Gauss 消去法

消元过程的第 $k$ 步，对矩阵需做 $(n-k)^2$ 次乘法运算及 $(n - k)$ 次除法运算，对右端向量需作 $(n - k)$ 次乘法运算，所以消去过程总的乘除法运算工作量为

$$\sum_{k=1}^{n-1}(n-k)^2+\sum_{k=1}^{n-1}(n-k)+\sum_{k=1}^{n-1}(n-k)=\frac13n^3+\frac12n^2-\frac56n.$$

回代过程中，计算每个 $x_k$ 需作 $(n - k + 1)$ 次乘除法运算，其工作量为

$$\sum_{k=1}^n(n-k+1)=\frac12n(n+1).$$

算法运行条件

要使得顺序消去法可运行，必须使得前 $n - 1$ 个主元素 $a^{(k)}_{kk} (k = 1，2， \dots，n - 1)$ 均不为零，其充要条件是系数矩阵 $A$ 的前 $n-1$ 个顺序主子式不为零，即

$$\left.D_k=\left|\begin{array}{ccc}a_{11}^{(1)}&\cdots&a_{1k}^{(1)}\\\vdots&&\vdots\\a_{k1}^{(1)}&\cdots&a_{kk}^{(1)}\end{array}\right.\right|\neq0,\quad k=1,2,\cdots,n-1.$$

为什么运行条件是主元素不为零充要是系数主子式不为零，这个其实也可以用LU分解的思路来思考

当主元素 $a^{(k)}_{kk}$ 均不为零时，我们可以考虑使用LU分解的思路来证明对应系数矩阵的顺序主子式不为零。

假设我们有一个 $n \times n$ 的方阵 $A$，其顺序主子式都不为零，即 $|A_k| \neq 0$，其中 $A_k$ 表示 $A$ 的第 $k$ 阶顺序主子式。

我们知道LU分解可以写为 $A = LU$，其中 $L$ 是单位下三角矩阵，$U$ 是上三角矩阵。我们可以将 $A$ 分解为：

$$A = \begin{bmatrix} a_{11} & 0 & \cdots & 0 \\ a_{21} & a_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{bmatrix} = \begin{bmatrix} 1 & 0 & \cdots & 0 \\ l_{21} & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ l_{n1} & l_{n2} & \cdots & 1 \end{bmatrix} \begin{bmatrix} u_{11} & u_{12} & \cdots & u_{1n} \\ 0 & u_{22} & \cdots & u_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & u_{nn} \end{bmatrix}$$

我们可以看到，对于LU分解中的 $U$ 矩阵，其对角线上的元素 $u_{kk}$ 对应于原矩阵 $A$ 的主元素 $a^{(k)}_{kk}$。由于主元素均不为零，因此LU分解中的 $U$ 矩阵的对角线上的元素 $u_{kk}$ 也均不为零。

现在我们来看LU分解中的 $L$ 矩阵。由LU分解的性质可知，$L$ 的对角线上的元素均为1，而且 $L$ 的第 $k$ 行第 $j$ 列的元素 $l_{kj}$ 是原矩阵 $A$ 的第 $j$ 阶顺序主子式与 $U$ 的第 $k$ 阶顺序主子式的商。由于 $A$ 的顺序主子式都不为零，因此 $L$ 的元素 $l_{kj}$ 也不为零。

一个问题

前面我们提到，除法运算中，要尽量避免除数的绝对值远小于被除数的绝对值

如果对角线元素 $a_{ii}$ 很小（导致 $|m_{ik}|$ 很大），就会导致很大的误差

列主元素Gauss 消去法

第 $k$ 次消元之前，通过行变换将 $a^{(k)}_{ik} (i =k，k + 1,\dots ，n)$ 中绝对值最大的元素交换到第 $k$ 行的主对角线位置上，然后进行消元计算，即：

其他过程同顺序法，算法要求各个列的主元素不为零

对于列主元素交换，相当于多乘了一个初等行变换矩阵，所以最终还是可以化为一个LU分解的形式

比如乘以

$$P = \begin{bmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix}$$

交换一二行

LU分解及其与Gauss消去的关系

考试不考，仅供拓展思路

对于 $Ax=b$ 将 $A$ 分解为 $A=LU$，其中 $L$ 是下三角矩阵， $U$ 是上三角矩阵，则可将原方程化为：

$$Ly=b,\,\,Ux=y$$

先由 $Ly=b$ 解出 $y$，再由 $Ux=y$ 解出 $x$

而Gauss消去也可以看作一个LU分解的过程，即：

$$L^{(n-1)}L^{(n-2)}\cdots L^1 \cdot A=A^{(n)}$$

其中 $L^1$ 对应于将第一列非主对角线元素都化为零的过程，即

$$L^1=\begin{bmatrix} 1 & 0 & 0 & \cdots &0\\ l_1 & 1 & 0 & \cdots &0\\ l_2 & 0 & 1 & \cdots &0\\ \vdots & \vdots & \vdots &&0\\ l_{n-1} & 0 & 0 & \cdots &1\\ \end{bmatrix}$$

其中 $l_i$ 对应 $-m_{i1}$ 。以此类推

不妨设 $L=(L^{(n-1)}L^{(n-2)}\cdots L^1)^{-1},U=A^{(n)}$，显然，$L$ 是下三角矩阵，$A$ 为上三角矩阵，即为 $A=LU$

矩阵条件数与病态方程组

设$A$、$\Delta A\in\mathbb{R}^n\times n$，$b$、$\Delta b\in\mathbb{R}^n$，$A$非奇异，$b\neq0$，$x$ 是方程组 $Ax=b$ 的解向量。若 $\|\Delta A\|<\frac1{\|{A}^{-1}\|}$ ，则有

• 方程组

  $$(A+\Delta A)(x+\Delta x)=b+\Delta b$$

  有唯一解 $x+\Delta x$

• 下列估计式成立：

  $$\frac{\parallel\Delta x\parallel}{\parallel x\parallel}\leqslant\frac{\parallel A\parallel\parallel A^{-1}\parallel}{1-\parallel A^{-1}\parallel\parallel\Delta A\parallel}(\frac{\parallel\Delta A\parallel}{\parallel A\parallel}+\frac{\parallel\Delta b\parallel}{\parallel b\parallel})$$

对于第一点，我们需要证明 $A+\Delta A$ 非奇异，证明关键其实是该定理：

设矩阵 $A \in R^{n\times n}$ 的某种范数 $\left\|A\right\|<1$ ，则 $I\pm A$ 为非奇异矩阵，并且当该种范数为算子范数时，下式成立

$$\left\|(I\pm A)^{-1}\right\|\leqslant\frac1{1-\|A\|}$$

由 $A+\Delta A=A(I+A^{-1}\Delta A)$，其中 $A$ 非奇异，即证 $I+A^{-1}\Delta A$ 非奇异，由上定理，即证 $\|A^{-1}\Delta A\|<1$

而 $\parallel A^{-1}\Delta A\parallel\leqslant\parallel A^{-1}\parallel\parallel\Delta A\parallel<1$，故得证

对于第二点，对原方程组进行展开，并代入 $Ax=b$

$$Ax+A\cdot \Delta x+\Delta A \cdot x+\Delta A\cdot \Delta x=b+\Delta b\\ A\cdot \Delta x+\Delta A \cdot x+\Delta A\cdot \Delta x=\Delta b\\ \Delta{x}={A}^{-1}\Delta{b}-{A}^{-1}\Delta{A}{x}-{A}^{-1}\Delta{A}\Delta{x}$$

两边同时取范数，由矩阵范数相容性，即 $\|b\|=\|Ax\|\leqslant\|A\|\|x\|$

$$\begin{aligned} \|\Delta x\| &= \|A^{-1}\Delta b\|+\|A^{-1}\Delta A\cdot x\|+\|A^{-1}\Delta{A}\Delta{x}\|\\ &\leq \parallel A^{-1}\parallel\parallel\Delta b\parallel+\parallel A^{-1}\parallel\parallel\Delta A\parallel\parallel x\parallel+\parallel A^{-1}\parallel\parallel\Delta A\parallel\parallel\Delta x\parallel \end{aligned}$$

即

$$(1-\parallel A^{-1}\parallel\parallel\Delta A\parallel)\frac{\parallel\Delta x\parallel}{\parallel x\parallel}\leqslant\parallel A\parallel\parallel A^{-1}\parallel\frac{\parallel\Delta b\parallel}{\parallel b\parallel}+\parallel A^{-1}\parallel\parallel\Delta A\parallel$$

由第一点 $\parallel A^{-1}\parallel\parallel\Delta A\parallel<1$，上式得证

一个细节就是这里用到的都是算子范数，使得 矩阵范数与向量范数相容，从而满足范数相容性

对于上面定理第二点的式子，我们可以改写成：

$$\frac{\parallel\Delta x\parallel}{\parallel x\parallel}\leqslant\frac{\parallel A\parallel\parallel A^{-1}\parallel}{1-\parallel A\parallel\parallel A^{-1}\parallel\frac{\parallel\Delta A\parallel}{\parallel A\parallel}}(\frac{\parallel\Delta A\parallel}{\parallel A\parallel}+\frac{\parallel\Delta b\parallel}{\parallel b\parallel})$$

那么 $\parallel A\parallel\parallel A^{-1}\|$ 越小，$A,b$ 的误差对解向量的相对误差影响就越小，反之则大。

条件数与病态

关于病态方程组求解

对于病态方程组的求解：

• 用更高精度的运算
• 平衡方法：$A$ 元素数量级差别很大时进行行平衡或列平衡
• 残差矫正法：病态但不特别严重时的一种缓解病态方法

迭代法

迭代法是目前求解大规模稀疏线性方程组的主要方法之一.

迭代法的一般形式及其收敛性

对于线性方程组，将系数矩阵分解为两个矩阵的差，即 $A=N-P$ ，代入 $Ax=b$，即：

$$Nx=Px+b\\x=N^{-1}Px+N^{-1}b$$

记 $G=N^{-1}P,\quad d=N^{-1}b$ 代入，得 $x=Gx+d$ ，取 $x^{(0)}\in\mathbb{R}^{n\times n}$ 作为初始近似解，使用递推式

$$x^{(k+1)}=Gx^{(k)}+d\quad(k=0,1,\cdots)$$

产生一个向量序列，当 $k$ 足够大时， $x^k$ 可以看作原方程的近似解

序列收敛性

只有在迭代法收敛的情况下，用它所产生的向量序列 $\{x^{(k)}\}$ 中的向量作为方程组的近似解才有意义，而且，k 越大，$x^{(k)}$ 作为方程组的解就越精确。关于向量序列收敛定义如下：

设有向量序列

$$x^{(k)}=\left({x_1}^{(k)},{x_2}^{(k)},\cdots,{x_n}^{(k)}\right)^\top\quad(k=0,1,\cdots)$$

如果存在常向量$x^*=(x_1^*,x_2^*,\cdots,x_n^*)^\top$，使得

$$\lim\limits_{k\to\infty}x_i^{(k)}=x_i^*\quad(i=1,2,\cdots,n)$$

则称向量序列$\{x^(k)\}$收敛于常向量$x^*$，记为

$$\lim\limits_{k\to\infty}x^{(k)}=x^*$$

相关定理

设有向量序列$\{x^{(k)}\}$和常向量 $x^*$，如果对某种范数，有

$$\lim\limits_{k\to\infty}\left\|x^{(k)}-x^*\right\|=0$$

则必有

$$\lim\limits_{k\to\infty}x^{(k)}=x^*$$

由向量范数等价性，原条件即：

$$\lim_{k\to\infty}\max_{1\leqslant i\leqslant n}\left|x_i^{(k)}-x_i^*\right|=0$$

一个后面用到的等式

$$x^{(k)}-x^*=\sum_{i=k}^\infty\left(x^{(i)}-x^{(i+1)}\right)$$

要证明这个公式，我们从Jacobi迭代法的定义出发。Jacobi迭代法用于求解线性方程组 $Ax = b$，迭代公式为：

$$x^{(k+1)} = D^{-1}(b - (L + U)x^{(k)})$$

其中 $A = D + L + U$，$D$ 是 $A$ 的对角矩阵，$L$ 是 $A$ 的严格下三角矩阵，$U$ 是 $A$ 的严格上三角矩阵。

我们要证明：$x^{(k)} - x^* = \sum_{i=k}^\infty (x^{(i)} - x^{(i+1)})$

首先，我们知道对于任何迭代步 $k$，有以下关系：

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

利用这一点，我们可以递推：

$$x^{(k+1)} - x^* = (x^{(k+1)} - x^{(k+2)}) + (x^{(k+2)} - x^*)$$

这样继续下去，我们得到：

$$x^{(k)} - x^* = (x^{(k)} - x^{(k+1)}) + (x^{(k+1)} - x^{(k+2)}) + (x^{(k+2)} - x^{(k+3)}) + \cdots$$

在极限的意义下：

$$x^{(k)} - x^* = \sum_{i=k}^\infty (x^{(i)} - x^{(i+1)})$$

这个公式成立的关键在于迭代过程的收敛性，即 $x^{(i)} \to x^*$ 当 $i \to \infty$。因此，尾项 $(x^{(i)} - x^*)$ 在极限下趋于零。

谱半径与收敛性

设 $n\times n$ 矩阵 $G$ 的特征值是 $\lambda_1,\lambda_2,\cdots,\lambda_n$ ,称

$$\rho(G)=\max_{1\leqslant i\leqslant n}|\lambda_i|$$

为矩阵 $G$ 的谱半径.

谱半径与收敛性的定理

• 定理：任意一个矩阵的谱半径与其范数是有关系的，即任一矩阵$A$的谱半径不大于$A$的任一与某一向量范数相容的矩阵范数，即

$$\rho(A)\leqslant\|A\|.$$

按照谱半径的定义 $\rho(A)=\max_{1\leq i\leq n}|\lambda_i|$， 其中，$\lambda_i(i=1,2,\cdots,n)$为 $n$ 阶方阵 $A$ 的特征值

设 $x$ 为对应于 $\lambda$ 的 $A$ 的特征向量，则有 $\lambda x=Ax$，两边取范数可得

$$|\lambda|\|x\|\leqslant\|A\|\|x\|.$$

因为 $x$ 为非零向量，$\|x\|\neq0$，故有$|\lambda|\leqslant\|A\|.$

上述这个不等式对$A$的任何特征值均成立，即可得 $\rho(A)\leqslant\|A\|.$

• 对上定理，利用矩阵的Jordan 标准型得：设 $A$ 是任意 $n$ 阶矩阵, 则 $A$ 的 $k$ 次幂 $A^k\to 0$ ( 当 $k \to +\infty$) 的充要条件为谱半径 $\rho(A)<1$

Jordan标准型是一种特殊形式的矩阵，通常用于简化矩阵的分析和计算特征值、特征向量等问题。对于一个给定的方阵，其Jordan标准型可以通过相似变换得到。以下是关于Jordan标准型的一些重要概念：

1. 定义：对于一个 $n \times n$ 的方阵 $A$，存在一个可逆矩阵 $P$，使得 $P^{-1}AP$ 是Jordan标准型。Jordan标准型是一个块对角矩阵，其中每个块称为Jordan块。
2. Jordan块：Jordan块是一个形式为 $\lambda I + N$ 的矩阵，其中 $\lambda$ 是矩阵的特征值，$I$ 是单位矩阵，$N$ 是上三角矩阵，其对角线元素为0或者都是1。Jordan块的大小取决于特征值的代数重数和几何重数。
3. 特征值和特征向量：Jordan标准型将矩阵的特征值和特征向量整合在一起，使得特征向量更易于计算和理解。每个Jordan块对应一个特征值，而特征向量则可以从Jordan块中推导出来。

• 定理：对任意的向量 $d$， 迭代法收敛的充分必要条件是 $\rho(G)<1$

对迭代形式 $x^{(k)}=Gx^{(k-1)}+d$ 和最终解 $x^*=Gx^*+d$ 做差并递推：

$$x^{(k)}-x^*=G\left(x^{(k-1)}-x^*\right)=\cdots=G^{k}\left(x^{(0)}-x^*\right)$$

若收敛，则 $x^{(k)}-x^*\to 0$ ，即 $G^{k}\left(x^{(0)}-x^*\right)\to 0$，即 $G^{k}\to 0,k\to \infty$，由上引理得 $\rho(G)<1$

同理可证充分性

收敛性判断

用迭代矩阵的谱半径判断迭代公式是否收敛往往不容易，可以使用如下定理

如果矩阵 $G$ 的某种范数$\|G\|<1$，则

• 方程组 $x=Gx+d$ 的解 $x^*$ 存在且唯一；

• 对于迭代公式 $x^{(k+1)}=Gx^{(k)}+d$ 有

  $$\lim\limits_{k\to\infty}x^{(k)}=x^*,\quad\forall x^{(0)}\in\mathbb{R}^n$$

并且下列两式成立

$$\left\|x^{(k)}-x^*\right\|\leqslant\frac{\|G\|^k}{1-\|G\|}\left\|x^{(1)}-x^{(0)}\right\| \tag{2.3.6}$$

$$\left\|x^{(k)}-x^*\right\|\leqslant\frac{\|G\|}{1-\|G\|}\left\|x^{(k)}-x^{(k-1)}\right\|\tag{2.3.7}$$

对于第一点，因 $\|G\| < 1$ , 根据定理算子范数的估计，可知矩阵 $I-G$非奇异， 其中 $I$ 是单位矩阵

对于第二点，由 $\|G\| < 1$ ，则 $\rho(G)<1$ ，即迭代法收敛，则 $\lim\limits_{k\to\infty}x^{(k)}=x^*$

则对 $x^{(k)}-x^*=\sum_{i=k}^\infty\left(x^{(i)}-x^{(i+1)}\right)$ 两边取范数，由$\|x+y\|\leq\|x\|+\|y\|$ 得：

$$\begin{aligned}\left\|x^{(k)}-x^*\right\|&\leqslant\sum_{i=k}^\infty\left\|x^{(i)}-x^{(i+1)}\right\|=\sum_{i=k}^\infty\left\|G^i\left(x^{(0)}-x^{(1)}\right)\right\|\\&\leqslant\sum_{i=k}^\infty\|G\|^i\left\|x^{(0)}-x^{(1)}\right\|\\&=\frac{\|G\|^k}{1-\|G\|}\left\|x^{(1)}-x^{(0)}\right\|,\end{aligned}$$

PS. $x^{(k)}-x^*=\sum_{i=k}^\infty\left(x^{(i)}-x^{(i+1)}\right)$ 这个等式上面证过

对于2.3.7

$$x^{(k)}-x^*=G\left(x^{(k-1)}-x^*\right)=G\left(x^{(k-1)}-x^{(k)}\right)+G\left(x^{(k)}-x^*\right).$$

则：

$$\begin{aligned}\left\|x^{(k)}-x^*\right\|&\leqslant\|G\|\left\|x^{(k-1)}-x^{(k)}\right\|+\|G\|\left\|x^{(k)}-x^*\right\|,\\&\left\|x^{(k)}-x^*\right\|\leqslant\frac{\|G\|}{1-\|G\|}\left\|x^{(k)}-x^{(k-1)}\right\|.\end{aligned}$$

Jacobi 迭代法

若系数矩阵 $A=(a_{ij})\in\mathbb{R}^{n\times n}$ 满足条件 $a_{ii}\neq0(i=1,2,\cdots,n)$，则将 $A$ 分解为 $A=D+L+U$

其中 $D$ 为对角矩阵，$L$ 为严格下三角矩阵，$U$ 为严格上三角矩阵

相当于 $D=diag(a_{ii})$ 且显然 $D^{-1}$ 存在

对于一般形式，取 $N=D,P=-(L+U)$ ，得迭代式：

$$x^{(k+1)}=-D^{-1}(L+U)x^{(k)}+D^{-1}b\quad(k=0,1,\cdots)$$

其中 $x^{(0)}\in\mathbb{R}^n$ 任取，迭代矩阵：

$$G_J=-D^{-1}(L+U)$$

此外，$d=D^{-1}b$

由 $D^{-1}=diag(a_{ii}^{-1})$，分量形式为：

$$x_i^{(k+1)}=(-\sum_{j=1}^na_{ij}x_j^{(k)}+b_i)/a_{ii}\quad(i=1,2,\ldots,n;k=0,1,\ldots)\\ x_i^{(k+1)}=-\sum_{j=1}^n\frac{a_{ij}}{a_{ii}}x_j^{(k)}+\frac{b_i}{a_{ii}}\quad(i=1,2,\ldots,n;k=0,1,\ldots)$$

• 由定理：对任意的向量 $d$， 迭代法收敛的充分必要条件是 $\rho(G)<1$，得： Jacobi 迭代法收敛的充分必要条件是 $\rho(G_J)<1$

• 由定理：矩阵 $G$ 的某种范数$\|G\|<1$，则方程组 $x=Gx+d$ 的解 $x^*$ 存在且唯一，得 ：如果 $\|G_{J}\|<1$, 则 Jacobi 迭代法收敛.

引理：

若矩阵 $A\in\mathbb{R}^{n\times n}$ 是主对角线按行 (或按列) 严格占优阵，则 $A$ 是非奇异矩阵.

主对角线按行 (或按列) 严格占优阵，即满足：

$$|a_{ii}|>\sum_{j=1\atop j\neq i}^n|a_{ij}|\:,\quad i=1,2,\cdots,n.$$

即主对角线元素大于该行非对角线元素绝对值之和

如果系数矩阵$A$是主对角线按行(或按列) 严格占优阵，则用Jacobi迭代法求解必收敛.

若 $A$ 是对称正定矩阵，则雅可比迭代收敛的充要条件是 $2D - A$ 也是对称正定矩阵.

Gauss-Seidel迭代法

在迭代法一般形式中，取$N=D+L,P=-U$，形成以下的迭代公式

$$x^{(k+1)}=-(D+L)^{-1}Ux^{(k)}+(D+L)^{-1}b\quad(k=0,1,\cdotp\cdotp\cdotp)$$

其中：

$$G_G=-(D+L)^{-1}U\\ d=(D+L)^{-1}b$$

使用时为了避免求 $D+L$ 的逆矩阵，可以通过这种方式改写：

$$\begin{aligned}&(D+L)x^{(k+1)}=-Ux^{(k)}+b\\&Dx^{(k+1)}=-Lx^{(k+1)}-Ux^{(k)}+b\\&x^{(k+1)}=-D^{-1}Lx^{(k+1)}-D^{-1}Ux^{(k)}+D^{-1}b\end{aligned}$$