工程数学 计算方法 第三章 线性方程组的数值解法

线性方程组的数值解法

迭代法

$$Ax=b⟹\overline{x}=B\overline{x}+\overline{f}⟹{\overline{x}}^{(k+1)}=B{\overline{x}}^{(k)}+\overline{f}$$

问题：
B=？（迭代格式？）
能使吗？（收敛条件？）
好使吗？（收敛速度？误差程度？）

逐次逼近。从一个或多个初始量出发，按照一定的计算格式获得方程组数值解的方法。简单实用。
收敛性是迭代法的前提。
适用于高次稀疏矩阵（常为200阶以上）。

Jacobi迭代法

将AX=b改写为X=BX+f的形式，建立雅可比方法的迭代格式：
$$Ax=b⟹x^{(k+1)}=B{x}^{(k)}+f$$
其中，B称为迭代矩阵。
思想与不动点迭代类似，使用前需要判别收敛性。

问题：怎么用（迭代格式）？什么时候能用（收敛条件）？计算到什么程度（迭代停止条件）？

迭代格式

第n式留xn，其他移到等号右边，写为 xn=求和 的形式。
$$\begin{gathered} \{\begin{array}{l}{a}_{11}{x}_1+a_{12}{x}_2+...+a_{1n}{x}_n=b_1\\ a_{21}{x}_1+a_{22}{x}_2+...+a_{2n}{x}_n=b_2\\ ...\\ a_{n1}{x}_1+a_{n2}{x}_2+...+a_{nn}{x}_n=b_n\end{array}\stackrel{a_{ii}\ne 0}{⟹}\{\begin{array}{l}{x}_1=\frac{1}{a_{11}}(-a_{12}{x}_2-...-a_{1n}{x}_n+b_1)\\ x_2=\frac{1}{a_{22}}(-a_{21}{x}_1-...-a_{2n}{x}_n+b_2)\\ ...\\ x_n=\frac{1}{a_{nn}}(-a_{n1}{x}_1-...-a_{n(n-1)}{x}_{n-1}+b_n)\end{array} \\ \begin{array}{rl}Ax=b& \iff (D+L+U)x=b\\ & \iff Dx=-(L+U)x+b\\ & \iff x=-D^{-1}(L+U)x+D^{-1}b\\ & \iff x=Bx+f\end{array} \end{gathered}$$
得Jacobi迭代公式：$X^{(k+1)}=-D^{-1}(L+U)X^{(k)}+D^{-1}b,k=0,1,...$，其中$B=-D^{-1}(L+U)$称为Jacobi迭代矩阵。

例：
$$\begin{gathered} \{\begin{array}{l}{a}_{11}{x}_1+a_{12}{x}_2+a_{13}{x}_3=b_1\\ a_{21}{x}_1+a_{22}{x}_2+a_{23}{x}_3=b_2\\ a_{31}{x}_1+a_{32}{x}_2+a_{33}{x}_3=b_3\end{array}\stackrel{a_{ii}\ne 0}{⟹}\{\begin{array}{l}{x}_1=\frac{1}{a_{11}}(-a_{12}{x}_2-a_{13}{x}_3+b_1)\\ x_2=\frac{1}{a_{22}}(-a_{21}{x}_1-a_{23}{x}_3+b_2)\\ x_n=\frac{1}{a_{33}}(-a_{31}{x}_1-a_{32}{x}_2+b_3)\end{array} \\ \stackrel{a_{ii}\ne 0}{⟹}\{\begin{array}{l}{x}_1^{(k+1)}=\frac{1}{a_{11}}(-a_{12}{x}_2^{(k)}-a_{13}{x}_3^{(k)}+b_1)\\ x_2^{(k+1)}=\frac{1}{a_{22}}(-a_{21}{x}_1^{(k)}-a_{23}{x}_3^{(k)}+b_2)\\ x_n^{(k+1)}=\frac{1}{a_{33}}(-a_{31}{x}_1^{(k)}-a_{32}{x}_2^{(k)}+b_3)\end{array} \end{gathered}$$

收敛条件⭐

误差：
$${\overrightarrow{e}}^{(k)}=B^k{\overrightarrow{e}}^{(0)}$$
得收敛的充分条件$\Vert B\Vert <1$。

向量：$\text{x}$是否收敛。向量收敛：分量各自收敛。

影响收敛性的因素：初值和方法本身。

Jacobi迭代法收敛三种判断方法：

1. A行对角占优；

2. 充分条件：$\Vert B\Vert <1$；

3. ⭐充要条件：$B^k\to 0⟺\rho (B)<1$。含参判断收敛性必须用特征值来算。

对角占优：矩阵$A=(a_{ij})\in {\text{R}}^{n\times n}$，若其满足$|a_{ij}|⩾{\sum }_{j=1,j\ne i}^n|a_{ij}|,i=1,2,...,n,$则称其对角占优。(行对角占优，即对角线上的数值比同行其他数绝对值之和大)

停止条件

一般精度要求为$||\text{X}-{\text{X}}^{(k)}|{|}_i⩽\epsilon$。因为$R^n$上范数等价，而且1-范数计算简单且更为严格，故实际计算中常取1-范数。且实际中无法得知精确值，故一般取两步之间的插值，即当
$$||{\text{X}}^{(k+1)}-{\text{X}}^{(k)}|{|}_1⩽\epsilon$$
时停止迭代。

Gauss-Seidel迭代法

对于Jacobi迭代法，每一步计算的时候都已经计算出新值，但靠后的分量仍然是用旧值计算的。优化：将式子中的值都使用最新值代入迭代，第i式中将$x_1$至$x_i-1$都使用第k+1次的迭代值代入计算。即：
$$\begin{array}{cc}{x}_i^{(k+1)}& =\frac{1}{a_{ii}}(-a_{i1}{x}_1^{(k+1)}-...-a_{i(i-1)}{x}_{i-1}^{(k+1)}-a_{i(i+1)}{x}_{i+1}^{(k)}-...-a_{in}{x}_n^{(k)}+b_i)\\ & =\frac{1}{a_{ii}}(b_i-{\sum }_{j=1}^{i-1}{a}_{ij}{x}_j^{(k+1)}-{\sum }_{j=i+1}^n{a}_{ij}{x}_j^{(k)})(i=1,2,...,n)\end{array}$$
矩阵形式：
$$\begin{array}{rl}{x}^{(k+1)}=-D^{-1}(L{x}^{(k+1)}+U{x}^{(k)})+D^{-1}b& \iff (D+L)x^{(k+1)}=-U{x}^{(k)}+b\\ & \iff x^{(k+1)}=-(D+L{)}^{-1}U{x}^{(k)}+(D+L{)}^{-1}b\\ & \iff x^{(k+1)}=B{x}^{(k)}+f\end{array}$$

其中$B=-(D+L{)}^{-1}U$称为Gauss-Seidel迭代矩阵。

与Jacobi迭代相比，这里每一次计算都是用的是最新计算得的数值。

可以原位运算。

收敛条件：

1. A行对角占优；
2. A为正定矩阵；
3. 充分条件：$\Vert B\Vert <1$；
4. 充要条件：$\rho (B)<1$。

判断正定矩阵常用条件：一切顺序主子式为正。i阶顺序主子式：取左上方i×i阶矩阵，这个矩阵对应的方程即为i阶顺序主子式。

精度判断与Jacobi迭代法相同。

SOR方法

逐次超松弛迭代法(Successive Over Relaxation Method)，课看作带参数ω的Gauss-Seidel迭代法。是对G-S迭代法的修正或加速。

算法构造

$$\begin{gathered} {\overline{x}}_i^{(k+1)}=\frac{1}{a_{ii}}(b_i-\sum _{j=1}^{i-1}{a}_{ij}{x}_j^{(k+1)}-\sum _{j=i+1}^n{a}_{ij}{x}_j^{(k)})(i=1,2,...,n) \\ \begin{array}{cc}{x}_i^{(k+1)}& =(1-\omega )x_i^{(k)}+\omega {\overline{x}}_i^{(k+1)}\\ & =x_i^{(k)}+\omega ({\overline{x}}_i^{(k+1)}-x_i^{(k)})\end{array} \end{gathered}$$

ω称为松弛因子。

由G-S迭代得修正前的值${\overline{x}}_i^{(k+1)}$，并通过这个式子来修正，得到迭代值$x_i^{(k+1)}$。

0<ω<1 :低松弛法

ω=1 :G-S迭代法

1<ω<2 :超松弛法

SOR迭代法的计算（不考这个计算）

（每个式子都是通式，代表一个方程组）

先移项：
$$x_i=\frac{1}{a_{ii}}(-a_{i1}{x}_1-...-a_{i(i-1)}{x}_{i-1}-a_{i(i+1)}{x}_{i+1}-...-a_{in}{x}_n+b_i)$$
写G-S迭代公式：
$$x_i^{(k+1)}=\frac{1}{a_{ii}}(-a_{i1}{x}_1^{(k+1)}-...-a_{i(i-1)}{x}_{i-1}^{(k+1)}-a_{i(i+1)}{x}_{i+1}^{(k)}-...-a_{in}{x}_n^{(k)}+b_i)$$
改写为SOR公式：
$$\begin{gathered} x_i^{(k+1)}=(1-\omega )x_i^{(k)}+\omega \frac{1}{a_{ii}}(-a_{i1}{x}_1^{(k+1)}-...-a_{i(i-1)}{x}_{i-1}^{(k+1)}-a_{i(i+1)}{x}_{i+1}^{(k)}-...-a_{in}{x}_n^{(k)}+b_i) \\ 即x_i^{(k+1)}=(1-\omega )x_i^{(k)}+\omega {\overline{x}}_i^{(k+1)} \end{gathered}$$
选取ω并带入初值。

直接法

直接法是理论上可以求得精确解的方法。

Guess消元法

1. 消元
2. 回代

矩阵：

1. 化上三角矩阵：前推过程
   主对角线上有0：换行/换列，$D_i\ne 0$时可换行
2. 自下而上解出自变量：回代过程

问题：化上三角矩阵时若有小主元会导致舍入误差的扩散，从而导致计算失败。
解决：全主元消去法和列主元消去法。全主元消去法更加稳定。

Gauss-Jordan消去法：主元变1，上下消为0。常用与求逆矩阵，不用于解方程。$[A|I]\Rightarrow [I|A^{-1}]$。

Doolittle分解法⭐

有时系数矩阵是确定的，而对于不同的问题只是b不同，而Guess消元法需要重新计算。

此外前推过程计算繁琐且易造成误差，而回代过程较为简单：能不能只回代不前推？或者增加回代的过程减少前推的过程：用回代过程代替前推过程。（回代过程的条件：A化为上三角矩阵/下三角矩阵。）

$A\overrightarrow{x}=\overrightarrow{b}$，A表示$x_i$之间的线性组合，表现了$x_i$间的相关性。若A定，分解成$LU\overrightarrow{X}=\overrightarrow{b}$，其中LU为三角阵。
$$LU\overrightarrow{X}=\overrightarrow{b}\Rightarrow \{\begin{array}{l}L\overrightarrow{y}=\overrightarrow{b}\\ U\overrightarrow{x}=\overrightarrow{y}\end{array}$$
分解然后二次迭代来避免前推过程。

可靠吗？

定理：A正定 $\iff$A的顺序主子式均不为0 $\Rightarrow$A的LU分解唯一

如何分解？

1. 乘积递推：（麻烦，运算量浪费大，不这么写）
   $$\begin{gathered} \{\begin{array}{l}{A}^{(n)}=L_1{L}_2...L_{n-2}{L}_{n-1}{A}^{(1)}\\ b^{(n)}=L_1{L}_2...L_{n-2}{L}_{n-1}{b}^{(1)}\end{array}⟹\{\begin{array}{l}A=L_1^{-1}{L}_2^{-1}...L_{n-2}^{-1}{L}_{n-1}^{-1}{A}^{(n)}\\ b=L_1^{-1}{L}_2^{-1}...L_{n-2}^{-1}{L}_{n-1}^{-1}{b}^{(n)}\end{array} \\ 则L\triangleq L_1^{-1}{L}_2^{-1}...L_{n-2}^{-1}{L}_{n-1}^{-1}为下三角矩阵，U=A^{(n)} \end{gathered}$$

2. 比较法

   根据A=LU导出L和U:
   $$\begin{gathered} \left[\begin{array}{cccc}{a}_{11}& \cdots & \cdots & a_{1n}\\ ⋮& & & ⋮\\ ⋮& & & ⋮\\ a_{n1}& \cdots & \cdots & a_{nn}\end{array}\right]=\left[\begin{array}{cccc}1& & & \\ l_{21}& 1& & \\ ⋮& \cdots & & \\ l_{n1}& \cdots & \cdots & 1\end{array}\right]\left[\begin{array}{cccc}{u}_{11}& \cdots & \cdots & u_{1n}\\ & & \cdots & ⋮\\ & & & ⋮\\ & & & u_{nn}\end{array}\right] \\ \{\begin{array}{ll}{u}_{1i}=a_{1i},& i=1,2,...,n,\\ l_{i1}=\frac{a_{i1}}{u_{11}},& i=2,3,...,n,\\ u_{ri}=a_{ri}-{\sum }_{k=1}^{r-1}{l}_{rk}{u}_{ki}& i=r,r+1,...,n,r=2,3,...,n,\\ l_{ri}=(a_{ir}-{\sum }_{k=1}^{r-1}{l}_{ik}{u}_{kr})/u_{rr}& i=r,r+1,...,n,r=2,3,...,n,\end{array} \end{gathered}$$
   u第一行，l第一列，u第二行，l第二列……（公式只是明确一点计算方法，实际计算之间按照这个顺序就能写出来了，写过程的时候也不需要写出分解的过程，直接得到UL就可以了）

   可以使用原位运算：$a_{ir}$只使用一次，用完后就不会在用到了，可把UL存在A的原地址中，节约空间。

扰动分析

实际求解中可能会有输入输出的微小扰动。不同方程组对扰动的敏感程度不同。分析对扰动的敏感程度：

定理：设线性方程组系数矩阵A及右端向量b分别带有微小扰动${\mathbf{\delta }}_{\mathbf{A}}$和${\mathbf{\delta }}_{\mathbf{b}}$，其导出的扰动解$\mathbf{X}+{\mathbf{\delta }}_{\mathbf{X}}$满足
$$(\mathbf{A}+{\mathbf{\delta }}_{\mathbf{A}})(\mathbf{X}+{\mathbf{\delta }}_{\mathbf{X}})=\mathbf{b}+{\mathbf{\delta }}_{\mathbf{b}}$$
则其解具有如下局部相对误差估计：
$$\frac{\Vert {\mathbf{\delta }}_{\mathbf{X}}\Vert }{\Vert \mathbf{X}\Vert }⩽\frac{\Vert {\mathbf{A}}^{-1}\Vert \Vert \mathbf{A}\Vert }{1-\Vert {\mathbf{A}}^{-1}\Vert \Vert {\mathbf{\delta }}_{\mathbf{A}}\Vert }\left(\frac{\Vert {\mathbf{\delta }}_{\mathbf{b}}\Vert }{\Vert \mathbf{b}\Vert }+\frac{\Vert {\mathbf{\delta }}_{\mathbf{A}}\Vert }{\Vert \mathbf{A}\Vert }\right)$$

当扰动${\mathbf{\delta }}_{\mathbf{A}}$充分小时，有近似关系：
$$\begin{gathered} \frac{\Vert {\mathbf{A}}^{-1}\Vert \Vert \mathbf{A}\Vert }{1-\Vert {\mathbf{A}}^{-1}\Vert \Vert {\mathbf{\delta }}_{\mathbf{A}}\Vert }\approx \Vert {\mathbf{A}}^{-1}\Vert \Vert \mathbf{A}\Vert \\ ⟹\frac{\Vert {\mathbf{\delta }}_{\mathbf{X}}\Vert }{\Vert \mathbf{X}+{\mathbf{\delta }}_{\mathbf{X}}\Vert }⩽\Vert {\mathbf{A}}^{-1}\Vert \Vert {\mathbf{\delta }}_{\mathbf{A}}\Vert =\Vert {\mathbf{A}}^{-1}\Vert \Vert \mathbf{A}\Vert \frac{\Vert {\mathbf{\delta }}_{\mathbf{A}}\Vert }{\Vert \mathbf{A}\Vert } \\ \therefore \frac{\Vert {\mathbf{\delta }}_{\mathbf{A}}\Vert }{\Vert \mathbf{A}\Vert }被放大了\Vert {\mathbf{A}}^{-1}\Vert \Vert \mathbf{A}\Vert 倍。 \end{gathered}$$
定义：设A为非奇异矩阵（满秩 / |A|=0），称数$Cond(A{)}_v=\Vert {\mathbf{A}}^{-1}{\Vert }_v\Vert \mathbf{A}{\Vert }_v(v=1,2或\infty )$为矩阵A的条件数。
C o n d 1 ( A ) ≜ ∥ A − 1 ∥ 1 ∥ A ∥ 1   C o n d ∞ ( A ) ≜ ∥ A − 1 ∥ ∞ ∥ A ∥ ∞   C o n d 2 ( A ) ≜ ∥ A − 1 ∥ 2 ∥ A ∥ 2 = max ⁡ 1 ⩽ i ⩽ n { λ i A T A } min ⁡ 1 ⩽ i ⩽ n { λ i A T A }   若 A 为 对 称 矩 阵 ， C o n d 2 ( A ) = max ⁡ ∣ λ A ∣ min ⁡ ∣ λ A ∣ Cond_1(A)\triangleq\Vert\boldsymbol{A}^{-1}\Vert_1\Vert\boldsymbol{A}\Vert_1\\\,\\ Cond_\infty(A)\triangleq\Vert\boldsymbol{A}^{-1}\Vert_\infty\Vert\boldsymbol{A}\Vert_\infty\\\,\\ Cond_2(A)\triangleq\Vert\boldsymbol{A}^{-1}\Vert_2\Vert\boldsymbol{A}\Vert_2= \sqrt{\frac{\underset{1\leqslant i\leqslant n}{\max}\{\lambda_i^{\boldsymbol{A}^T\boldsymbol{A}}\}} {\underset{1\leqslant i\leqslant n}{\min}\{\lambda_i^{\boldsymbol{A}^T\boldsymbol{A}}\}}}\\\,\\ 若A为对称矩阵， Cond_2(A)=\sqrt{\frac{\max\vert\lambda^A\vert}{\min\vert\lambda^A\vert}} Cond1​(A)≜∥A−1∥1​∥A∥1​Cond∞​(A)≜∥A−1∥∞​∥A∥∞​Cond2​(A)≜∥A−1∥2​∥A∥2​=1⩽i⩽nmin​{λiATA​}1⩽i⩽nmax​{λiATA​}​ ​若A为对称矩阵，Cond2​(A)=min∣λA∣max∣λA∣​ ​
当$Cond(A)>>1$时，视为方程组是病态的，即当输入有微小偏差时扰动解$\hat{\mathbf{X}}$存在较大误差。

病态与否为方程组的性质，病态矩阵没有任何补救方法，只能修改计算方法。（这一部分记这一句话就行了）

应试

迭代法的收敛性

收敛性完全取决于迭代矩阵的性质，与初值选取无关。
考试可能只考二阶，三阶不好求逆矩阵（线代忘光）。
迭代的求解计算量大，不考，只计算收敛性。记住迭代公式的形式，记矩阵形式即可，然后用$\rho (B)<1$ 来判断收敛性。
$\rho (B)$ ：$B$ 特征值绝对值的最大值。

迭代格式里要求逆矩阵：
二阶矩阵：主交换，副相反，除以行列式
三阶矩阵：

1. 初等变换法
2. 伴随式除以行列式
   伴随式：伴随阵中$a_{ij}$为原矩阵第$j$行和第$i$列划掉，剩下的部分求行列式的值，并乘上$(-1{)}^{i+j}$。例：$a_{22}=(-1{)}^4\left|\begin{array}{cc}{a}_{11}& a_{13}\\ a_{31}& a_{33}\end{array}\right|$ 。

求$\rho (B)$：求矩阵特征值：
解方程：$|\lambda I-B|=0$ 。若为三阶矩阵注意因式分解，先找一个因式才好解。

Doolittle分解法

考试考三阶的。
分解方法：u第一行，l第一列，u第二行，l第二列，u第三行，l第三列。

可能会考条件数的计算，和矩阵范数一起考。概念在上文扰动分析部分。