数值计算笔记-线性方程组迭代解法

三种经典迭代法python例子

1. 矩阵极限： { A ( k ) } \{A^{(k)}\} {A(k)}为矩阵的序列，且$A=(a_{ij})\in R^{n\times n}$，若有
$$\underset{k\to \infty }{\lim }{a}_{ij}^{(k)}=a_{ij},(i,j=1,\cdots ,n)$$
则 { A ( k ) } \{A^{(k)}\} {A(k)}收敛于A，记为${\lim }_{k\to \infty }{A}^{(k)}=A$

算子范数等价性：$\exists c_1,c_2>0$对$\forall A\in R^{n\times n}$，和任意算子范数$||\cdot |{|}_i$，有
$$c_1||A|{|}_{\infty }\le ||A|{|}_i\le c_2||A|{|}_{\infty }$$
范数与极限：其中$||\cdot ||$表示任意一种算子范数
$$\underset{k\to \infty }{\lim }{A}^{(k)}=A\iff \underset{k\to \infty }{\lim }||A^{(k)}-A||=0$$
极限与线性方程：${\lim }_{k\to \infty }{A}^{(k)}=A$充要条件：$\forall \vec{x}$有${\lim }_{k\to \infty }{A}^{(k)}\vec{x}=A\vec{x}$

2. 定义谱半径（ = 最大特征值）：$A\in R^{n\times n}$特征值为${\lambda }_i,(i=1,\cdots ,n)$，则 ρ ( A ) = max ⁡ 1 < i ≤ n { λ i } \rho(A) = \max_{1 < i \leq n} \{\lambda_i\} ρ(A)=max1<i≤n​{λi​}成为谱半径

谱半径与算子范数：

（1）$\rho (A)\le ||A||$，即谱半径不大于任意算子范数

（2）$||A|{|}_2=\rho (A)$

谱半径与极限：$B=(b_{ij})\in R^{n\times n}$，则${\lim }_{k\to \infty }{B}^{n\times n}=0\iff \rho (B)<1$

3. 迭代法基本形式：原方程$A\vec{x}=\vec{b}$，迭代形式${\vec{x}}^{(k+1)}=B{\vec{x}}^{(k)}+\vec{f}$

算法：1阶定长迭代法

输入：初值${\vec{x}}_0$，B，$\vec{f}$；输出：$\vec{x}$
$$\begin{array}{rl}& \vec{x}={\vec{x}}_0\\ & while()\\ & \vec{x}=B\vec{x}+\vec{f}\\ & End\end{array}$$
4. 迭代法全局收敛充要条件：若迭代法${\vec{x}}_{k+1}=B{\vec{x}}_k+\vec{f}$，中I - B非奇异。对$\forall {\vec{x}}_0$迭代法得到的序列 { x ⃗ k } \{\vec{x}_k\} {x k​}收敛$\iff \rho (B)<1$。在这种条件下 { x ⃗ k } \{\vec{x}_k\} {x k​}的极限${\vec{x}}^{\ast }$为$\vec{x}=B\vec{x}+\vec{f}$唯一解

收敛条件：迭代形式${\vec{x}}_{k+1}=B{\vec{x}}_k+\vec{f}$，若$\forall ||B||=q<1$，则

（1）全局收敛

（2）$||{\vec{x}}_k-{\vec{x}}^{\ast }||\le q^k||{\vec{x}}_0-{\vec{x}}^{\ast }||$

（3）$||{\vec{x}}_k-{\vec{x}}^{\ast }||\le \frac{q}{1-q}||{\vec{x}}_k-{\vec{x}}_{k+1}||$

（4）$||{\vec{x}}_k-{\vec{x}}^{\ast }||\le \frac{q^k}{1-q}||{\vec{x}}_1-{\vec{x}}_0||$

定义收敛阶：解序列为${\vec{x}}_k,{\vec{x}}_k\in R^n$，收敛到向量${\vec{x}}^{\ast }$，若
$$\begin{gathered} \underset{k\to \infty }{\lim }\frac{{\vec{e}}_{k+1}}{{\vec{e}}_k}=c\ne 0 \\ {\vec{e}}_k={\vec{x}}_k-{\vec{x}}^{\ast } \end{gathered}$$
则迭代过程是p阶收敛的。1阶定常迭代法收敛速度：${\lim }_{k\to \infty }\frac{||{\vec{e}}_{k+1}||}{||{\vec{e}}_k||}=\rho (B)$

定义收敛速度：$R=-lo{g}_{10}\rho (B)$

5. 雅可比迭代法

思路：第k次迭代得到的n维向量${\vec{x}}_k=(x_1^{(k)},\cdots ,x_n^{(k)})$。以$3\times 3$维矩阵A为例
$$\{\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}\Rightarrow \{\begin{array}{l}{x}_1^{}=-\frac{1}{a_{11}}(a_{12}{x}_2+a_{13}{x}_3)+\frac{b_1}{a_{11}}\\ x_2=-\frac{1}{a_{22}}(a_{21}{x}_1+a_{13}{x}_2)+\frac{b_2}{a_{22}}\\ x_3=-\frac{1}{a_{33}}(a_{31}{x}_1+a_{32}{x}_2)+\frac{b_3}{a_{33}}\end{array}$$
迭代公式：
$$\begin{array}{rl}& \{\begin{array}{l}A\vec{x}=\vec{b}\\ A=D-L-U\end{array}\\ & \Rightarrow {\vec{x}}_{k+1}=D^{-1}(L+U){\vec{x}}_k+D^{-1}\vec{b}\\ & \Rightarrow \{\begin{array}{l}{\vec{x}}_{k+1}=B{\vec{x}}_k+\vec{f}\\ B=D^{-1}(L+U)\\ \vec{f}=D^{-1}\vec{b}\end{array}\end{array}$$
其中D为A的对角线元素组成的对角矩阵；L、U为A下三角、上三角元素组成的元素，其对角元素为0

算法：输入：初值$\vec{x}$，A，$\vec{b}$；输出：$\vec{x}$
$$\begin{array}{rl}& while(不满足判停条件)\\ & \vec{y}=\vec{x}\\ & fori=1\to n\\ & x_i=(b_i-\sum _{j=1}^{i-1}{a}_{ij}{y}_j-\sum _{j=i+1}^n{a}_{ij}{y}_j)/a_{ii}\\ & End\\ & End\end{array}$$
6.高斯-塞尔德迭代法（G-S迭代法）

思路：
$$\{\begin{array}{l}{x}_1^{(k+1)}=-\frac{1}{a_{11}}(+a_{12}{x}_2^{(k)}+a_{13}{x}_3^{(k)})+\frac{b_1}{a_{11}}\\ x_2^{(k+1)}=-\frac{1}{a_{22}}(a_{21}{x}_1^{(k+1)}+a_{23}{x}_3^{(k)})+\frac{b_2}{a_{22}}\\ x_3^{(k+1)}=-\frac{1}{a_{33}}(a_{31}{x}_1^{(k+1)}+a_{32}{x}_2^{(k+1)})+\frac{b_3}{a_{3}3}\end{array}$$
迭代公式
$$\begin{array}{rl}& \{\begin{array}{l}A\vec{x}=\vec{b}\\ A=D-L-U\end{array}\\ & \Rightarrow {\vec{x}}_{k+1}=D^{-1}(L{\vec{x}}_{k+1}+U{\vec{x}}_k)+D^{-1}\vec{b}\\ & \Rightarrow {\vec{x}}_{k+1}=(D-L{)}^{-1}U{\vec{x}}_k+(D-L{)}^{-1}\vec{b}\\ & \Rightarrow \{\begin{array}{l}{\vec{x}}_{k+1}=B{\vec{x}}_k+\vec{f}\\ B=(D-L{)}^{-1}U\\ \vec{f}=(D-L{)}^{-1}\vec{b}\end{array}\end{array}$$
算法：输入：初值$\vec{x}$，A，$\vec{b}$；输出：近似解$\vec{x}$
$$\begin{array}{rl}& while(不满足判停条件)\\ & fori=1\to n\\ & x_i=(b_i-\sum _{j=1}^{i-1}{a}_{ij}{x}_j-\sum _{j=i+1}^n{a}_{ij}{x}_j)/a_{ii}\\ & End\\ & End\end{array}$$
7. 逐次超松弛迭代法（SOR迭代法）

思路：由G-S迭代法得到${\vec{x}}_{k+1}^{\prime }$。则SOR迭代得到的向量为${\vec{x}}_{k+1}=(1-w){\vec{x}}_k+w{\vec{x}}_{k+1}^{\prime }$

迭代公式：
$$\begin{array}{rl}& \{\begin{array}{l}A\vec{x}=\vec{b}\\ A=D-L-U\end{array}\\ & \Rightarrow {\vec{x}}_{k+1}=(1-w){\vec{x}}_k+w[D^{-1}(L{\vec{x}}_{k+1}+U{\vec{x}}_k)]+D^{-1}\vec{b}\\ & \Rightarrow {\vec{x}}_{k+1}=(D-wL{)}^{-1}((1-w)D+wU){\vec{x}}_k+(D-sL{)}^{-1}w\vec{b}\\ & \Rightarrow \{\begin{array}{l}{\vec{x}}_{k+1}=B{\vec{x}}_k+\vec{f}\\ B=(D-wL{)}^{-1}((1-w)D+wU)\\ \vec{f}=(D-sL{)}^{-1}w\vec{b}\end{array}\end{array}$$
算法：输入：初值$\vec{x}$，A，$\vec{b}$，参数w；输出：近似解$\vec{x}$
$$\begin{array}{rl}& while(不满足判停条件)\\ & fori=1\to n\\ & x_i=(1-w)x_i+w(b_i-\sum _{j=1}^{i-1}{a}_{ij}{x}_j-\sum _{j=i+1}^n{a}_{ij}{x}_j)/a_{ii}\\ & End\\ & End\end{array}$$

8. 雅可比迭代收敛的充要条件：已知对角元素$a_{ii}>0(i=1,\cdots ,n)$。$A\vec{x}=\vec{b}$迭代法收敛$\iff$A和2D-A正定，其中D为A对角线元素组成的矩阵

推论：设B为雅可比迭代中矩阵。若$||B|{|}_{\infty }<1$或者$||B|{|}_1<1$，则G-S迭代收敛

9. 定义可约矩阵：设$A=(a_{ij})\in R^{n\times n}$，若有排列矩阵P（一系列初等矩阵的乘积）使
$$P^{T}AP=\left[\begin{array}{cc}{A}_{11}& A_{12}\\ O& A_{22}\end{array}\right]$$
其中$A_{11},A_{22}$为阶数不等于0的方阵，则A为可约矩阵

对角占优性质：若A为严格对角占优矩阵，或不可约的弱对角占优矩阵，则A非奇异

10. G-S/SOR迭代法收敛条件

（1）若矩阵A严格对角占优，或不可约对角占优，则雅可比迭代法和G-S迭代法收敛

（2）条件同（1），松弛因子$0<w\le 1$的SOR迭代法收敛

正定矩阵收敛性质

（1）若A对称正定，则G-S迭代法和$0<w<2$的SOR迭代法收敛

（2）若SOR迭代法收敛，则$w\in (0,2)$

11. 最速下降法（= 梯度下降法）

推导：

（1）用变分原理将$A\vec{x}=\vec{b}$变为求$\varphi (\vec{x})=\frac{1}{2}{\vec{x}}^{T}A\vec{x}-{\vec{b}}^T\vec{x}$最小值点

（2）迭代方式：${\vec{x}}_{k+1}={\vec{x}}_k+{\alpha }_k{\vec{p}}_k$，其中${\alpha }_k$搜索步长，${\vec{p}}_k$搜索方向

（3）确定方向：下降速度最开的方向为负梯度方向，即$-\nabla \varphi (\vec{x})=-(\frac{\partial \varphi }{\partial x_1},\cdots ,\frac{\partial \varphi }{\partial x_n}{)}^T$

（4）确定步长
$$\begin{array}{rl}{\alpha }_k& =arg\underset{\alpha }{\min }\frac{1}{2}({\vec{x}}_k+\alpha {\vec{p}}_k{)}^{T}A({\vec{x}}_k+\alpha {\vec{p}}_k)-{\vec{b}}^T({\vec{x}}_k+\alpha {\vec{p}}_k)\\ & =arg\underset{\alpha }{\min }\frac{1}{2}{\alpha }^2{\vec{p}}_k^{T}A{\vec{p}}_k-\alpha {\vec{r}}_k{\vec{p}}_k+\varphi ({\vec{x}}_k)({\vec{r}}_k=\vec{b}-A\vec{x})\\ & \Rightarrow {\alpha }_k=\frac{{\vec{r}}_k^T{\vec{p}}_k}{{\vec{p}}_{k}A{\vec{p}}_k}\end{array}$$
算法：输入：初值$\vec{x}$，A，$\vec{b}$；输出：近似值解$\vec{x}$
$$\begin{array}{rl}& \vec{r}=\vec{b}-A\vec{x}\\ & while不满足判停条件\\ & \alpha ={\vec{r}}^T\vec{r}/({\vec{r}}^{T}A\vec{r})\\ & \vec{x}=\vec{x}+\alpha \vec{r}\\ & \vec{r}=\vec{r}-\alpha A\vec{r}\\ & End\end{array}$$

12. 共轭梯度法

思路：第一步与最速下降法一样，之后每次迭代的搜索方向由2个向量线性组合生成，即${\vec{x}}_{k+1}={\vec{x}}_k+\alpha (\xi {\vec{r}}_k+\eta {\vec{p}}_{k-1})$

推导搜索方向：设$f(\xi ,\eta )=\varphi ({\vec{x}}_k+\xi {\vec{r}}_k+\eta {\vec{p}}_{k-1})$

令$\frac{\partial f}{\partial \xi }=\frac{\partial f}{\partial \eta }=0$，得到
$$\begin{gathered} \{\begin{array}{l}\xi {\vec{r}}_k^{T}A{\vec{x}}_r+\eta {\vec{r}}_k^{T}A{\vec{p}}_{k-1}={\vec{r}}^T\vec{r}\\ \xi {\vec{r}}_k^{T}A{\vec{p}}_{k-1}+\eta {\vec{p}}_{k-1}^{T}A{\vec{p}}_{k-1}=0\end{array} \\ \Rightarrow {\vec{p}}_k=\frac{1}{\xi }({\vec{x}}_{k+1}-{\vec{x}}_k)={\vec{r}}_k+\frac{\eta }{\xi }{\vec{p}}_{k-1}={\vec{r}}_k+\beta {\vec{p}}_{k-1} \\ {\beta }_{k-1}=\frac{\eta }{\xi }=\frac{{\vec{r}}_k^{T}A{\vec{p}}_{k-1}}{{\vec{p}}_{k-1}^{T}A{\vec{p}}_{k-1}} \end{gathered}$$
步长推导不变，最终得到${\vec{x}}_{k+1}={\vec{x}}_k+{\alpha }_k{\vec{p}}_k$

算法：输入：初值${\vec{x}}_0$，A，$\vec{b}$；输出：${\vec{x}}_k$
$$\begin{array}{rl}& {\vec{r}}_0=\vec{b}-A{\vec{x}}_0\\ & {\vec{p}}_0={\vec{r}}_0\\ & k=0\\ & while不满足判停条件\\ & {\alpha }_k={\vec{r}}_k^T{\vec{p}}_k/({\vec{p}}_k^{T}A{\vec{r}}_k)\\ & {\vec{x}}_{k+1}={\vec{x}}_k+{\alpha }_k{\vec{p}}_k\\ & {\vec{r}}_{k+1}={\vec{r}}_k-{\alpha }_{k}A{\vec{p}}_k//下面计算下一次的步长和方向\\ & {\beta }_k=-{\vec{r}}_{k+1}^{T}A{\vec{p}}_k/({\vec{p}}_k^{T}A{\vec{p}}_k)\\ & {\vec{p}}_{k+1}={\vec{r}}_{k+1}+{\beta }_k{\vec{p}}_k\\ & k=k+1\\ & End\end{array}$$