数值分析：线性方程组迭代法

迭代法思想

迭代改善算法是求解病态矩阵的一种比较理想的方法。

对于给定的方程组

{a11⋅x1+a12⋅x2+⋯+a1n⋅xn=b1a21⋅x1+a22⋅x2+⋯+a2n⋅xn=b2⋯an1⋅x1+an2⋅x2+⋯+ann⋅xn=bn\begin{cases} 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 \\ \cdots \\ a_{n1} ⋅ x_1 + a_{n2} ⋅ x_2 + \cdots + a_{nn} ⋅ x_n = b_n \end{cases}

将其改写为x=B0x+fx=B_0x+f的形式

{x1=−(a12⋅x2+⋯+a1n⋅xn)a11+b1a11x2=−(a21⋅x1+⋯+a2n⋅xn)a22+b2a22⋯xn=−(an1⋅x1+an2⋅x2+⋯ )ann+bnann\begin{cases} x_1 = - \dfrac{ (a_{12} ⋅ x_2 + \cdots + a_{1n} ⋅ x_n) }{ a_{11} } + \dfrac{b_1}{a_{11}} \\ x_2 = - \dfrac{ (a_{21} ⋅ x_1 + \cdots + a_{2n} ⋅ x_n) }{ a_{22} } + \dfrac{b_2}{a_{22}} \\ \cdots \\ x_n = - \dfrac{ (a_{n1} ⋅ x_1 + a_{n2} ⋅ x_2 + \cdots) }{ a_{nn} } + \dfrac{b_n}{a_{nn}} \end{cases}

B0=[0−a12a11⋯−a1na11−a21a220⋯−a2na22⋮⋮⋱⋮−an1ann−an2ann⋯0]B_0= \left[\begin{array}{c} 0 & -\frac{a_{12}}{a_{11}} & \cdots & -\frac{a_{1n}}{a_{11}} \\ -\frac{a_{21}}{a_{22}} & 0 & \cdots & -\frac{a_{2n}}{a_{22}} \\ \vdots & \vdots & \ddots & \vdots \\ -\frac{a_{n1}}{a_{nn}} & -\frac{a_{n2}}{a_{nn}} & \cdots & 0 \end{array}\right]，f=[b1a11b2a22⋮bnann]f= \left[\begin{array}{c} \frac{b_1}{a_{11}} \\ \frac{b_2}{a_{22}} \\ \vdots \\ \frac{b_n}{a_{nn}} \end{array}\right]

任取初值x(0)=[x1(0),x2(0),⋯ ,xn(0)]Tx^{(0)} = \left[ x_1^{(0)}, x_2^{(0)}, \cdots, x_n^{(0)} \right]^T，带入上式右端即可得：

x(1)=[x1(1),x2(1),⋯ ,xn(1)]Tx^{(1)} = \left[ x_1^{(1)}, x_2^{(1)}, \cdots, x_n^{(1)} \right]^T

重复上步，即可得：

x(k)=[x1(k),x2(k),⋯ ,xn(k)]T,   k∈Nx^{(k)} = \left[ x_1^{(k)}, x_2^{(k)}, \cdots, x_n^{(k)} \right]^T, \ \ \ k∈N

若lim⁡k→∞x(k)=x∗\lim\limits_{k→∞}x^{(k)} = x^*则称此迭代收敛。该方法称为迭代法（或一阶定常迭代法）。

迭代格式

x(k+1)=Bxk+f,   k∈Nx^{(k+1)}=Bx^{k}+f, \ \ \ k∈N

误差向量

为研究x(k)x^{(k)}的收敛性，引进误差向量ε(k+1)=x(k+1)−x∗ε^{(k+1)} = x^{(k+1)} - x^*

ε(k)=Bε(k−1)=B2ε(k−2)=⋯=Bkε(0)ε^{(k)} = Bε^{(k-1)} = B^{2}ε^{(k-2)} = \cdots = B^kε^{(0)}

迭代法收敛其实就是ε(k)ε^{(k)}收敛到0。当迭代序列收敛时，迭代序列的收敛解就是原方程的精确解。

迭代矩阵

设Ax=b,   ∣A∣≠0Ax=b, \ \ \ |A|≠0，将AA分裂为A=M−NA=M-N，其中MM称为分裂矩阵。

(M−N)x=b(M-N)x=b
   ⇒   Mx=Nx+b{ \ \ \ ⇒ \ \ \ } Mx=Nx+b
   ⇒   x=M−1Nx+M−1b{ \ \ \ ⇒ \ \ \ } x=M^{-1}Nx+M^{-1}b

x=(M−1N)⏟x+(M−1b)⏟Bf\begin{matrix} x &=& \underbrace{(M^{-1}N)}x &+& \underbrace{(M^{-1}b)} \\ && B && f \end{matrix}

构造一阶定常迭代法

x(k+1)=Bxk+fx^{(k+1)}=Bx^{k}+f

A=M−N,  B=M−1NA=M-N ,\,\, B=M^{-1}N
  ⇒  B=M−1(M−A)\,\,⇒\,\, B = M^{-1}(M-A)
  ⇒  B=I−M−1A\,\,⇒\,\, B = I - M^{-1}A

B=I−M−1AB = I - M^{-1}A

同理，若将AA分解为A=M+NA=M+N则上式依然成立。

通过将AA分解为A=M−NA=M-N或A=M+NA=M+N，可构造出

{B=I−M−1Af=M−1b\begin{cases} B = I - M^{-1}A \\ f = M^{-1}b \end{cases}

迭代改善算法

①A=D−(L+U)②A=D+(L+U)\begin{matrix} ① & A = D - (L + U) \\ ② & A = D + (L + U) \end{matrix}

①和②均为正确构造方法，上课时，老师使用的是②。

雅可比迭代（Jacobi）

同步迭代算法

若aii≠0a_{ii}≠0，记①②{M=DN=L+U①② \begin{cases} M = D \\ N = L + U \end{cases}，则

迭代矩阵

B=I−M−1A={①I−D−1(D−L−U)=D−1(L+U)②I−D−1(D+L+U)=−D−1(L+U)B = I-M^{-1}A = \begin{cases} ① & I-D^{-1}(D-L-U) = D^{-1}(L+U) \\ ② & I-D^{-1}(D+L+U) = -D^{-1}(L+U) \end{cases}

计算公式

由①①得

Dx(k+1)=(L+U)x(k)+bDx^{(k+1)} = (L+U)x^{(k)} + b

记x(k)=[x1(k),xi(k),⋯ ,xn(k)]Tx^{(k)} = \left[ x_1^{(k)}, x_i^{(k)}, \cdots, x_n^{(k)} \right]^T

aiixi(k+1)=−∑j=1i−1aijxj(k)−∑j=i+1naijxj(k)+bia_{ii}x_i^{(k+1)} = -\sum\limits_{j=1}^{i-1} a_{ij} x_j^{(k)} -\sum\limits_{j=i+1}^{n} a_{ij} x_j^{(k)} + b_i

xi(k+1)=bi−∑j≠ij=1naijxj(k)aiix_i^{(k+1)} = \dfrac{ b_i - \sum\limits_{^{j=1}_{j≠i}}^{n} a_{ij}x_j^{(k)} }{ a_{ii} }

高斯-赛德尔迭代（Gauss-Seide）

异步迭代算法

①{M=D−LN=U① \begin{cases} M = D - L \\ N = U \end{cases}，②{M=D+LN=U② \begin{cases} M = D + L \\ N = U \end{cases}

迭代矩阵

B=I−M−1A={①I−M−1(M−U)=(D−L)−1U②I−M−1(M+U)=−(D+L)−1UB = I-M^{-1}A = \begin{cases} ① & I-M^{-1}(M-U) = (D-L)^{-1}U \\ ② & I-M^{-1}(M+U) = -(D+L)^{-1}U \end{cases}

计算公式

由①①得

Dx(k+1)=Dx(k)−(Lx(k+1)+Ux(k)−Dx(k)+b)Dx^{(k+1)} = Dx^{(k)} - (Lx^{(k+1)} + Ux^{(k)} - Dx^{(k)} + b)

xi(k+1)=bi−∑j=1i−1aijxj(k+1)−∑j=i+1naijxj(k)aiix_i^{(k+1)} = \dfrac{ b_i - \sum\limits_{j=1}^{i-1} a_{ij}x_j^{(k+1)} - \sum\limits_{j=i+1}^{n} a_{ij}x_j^{(k)} }{ a_{ii} }

例

分别用雅可比迭代法及高斯-赛德尔迭代法解下列方程组

{5x1+2x2+x3=−12−x1+4x2+2x3=202x1−3x2+10x3=3\begin{cases} 5x_1+2x_2+x_3=−12 \\ −x_1+4x_2+2x_3=20 \\ 2x_1−3x_2+10x_3=3 \end{cases}，取x(0)=[000]x^{(0)} = \left[\begin{array}{c} 0 \\ 0 \\ 0 \end{array}\right]

问

解

事前估计：k>ln⁡ε(1−q)∥x(1)−x(0)∥ln⁡qk > \dfrac{ \ln \dfrac{ε(1-q)}{\|x^{(1)}-x^{(0)}\|} }{ \ln q }

迭代构造：{A=D+L+UB=I−M−1A={−D−1(L+U)MJacobi=D−(D+L)−1UMGauss−Seide=D+Lf=M−1b\begin{cases} A = D + L + U \\ B = I - M^{-1}A = \begin{cases} -D^{-1}(L+U) & M_{Jacobi} = D \\ -(D+L)^{-1}U & M_{Gauss-Seide} = D + L \end{cases} \\ f = M^{-1}b \end{cases}

严格对角占优阵：∣aii∣>∑j≠ij=1n∣aij∣|a_{ii}| > \sum\limits_{^{j=1}_{j≠i}}^{n} |a_{ij}|

A=[521−1422−310]=L+D+U=[−12−3]+[5410]+[212]A = \left[\begin{array}{c} 5 & 2 & 1 \\ −1 & 4 & 2 \\ 2 & −3 & 10 \end{array}\right] = L + D + U = \left[\begin{array}{c} \\ −1 \\ 2 & −3 \end{array}\right] + \left[\begin{array}{c} 5 \\ & 4 \\ && 10 \end{array}\right] + \left[\begin{array}{c} & 2 & 1 \\ && 2 \\ & \end{array}\right]

Jacobi

{BJacobi=I−D−1A=[0−2/5−1/51/40−1/2−1/53/100]fJacobi=D−1b=[−1255310]T∥BJacobi∥∞=34<1\begin{cases} B_{Jacobi} = I - D^{-1}A = \left[\begin{array}{c} 0 & −2/5 & −1/5 \\ 1/4 & 0 & −1/2 \\ −1/5 & 3/10 & 0 \end{array}\right] \\ \\ f_{Jacobi} = D^{-1}b = \left[\begin{array}{c} -\dfrac{12}{5} & 5 & \dfrac{3}{10} \end{array}\right]^T \\ \\ \|B_{Jacobi}\|_{\infty} = \dfrac{3}{4} < 1 \end{cases}

x(1)=BJacobix(0)+fJacobi=[−1255310]Tx^{(1)} = B_{Jacobi}x^{(0)} + f_{Jacobi} = \left[\begin{array}{c} -\dfrac{12}{5} & 5 & \dfrac{3}{10} \end{array}\right]^T

∥x(1)−x(0)∥∞=5\|x^{(1)}-x^{(0)}\|_{\infty} = 5

k>ln⁡ε(1−q)∥x(1)−x(0)∥ln⁡q=ln⁡10−4(1−34)5ln⁡34≈42.43k > \dfrac{ \ln \dfrac{ε(1-q)}{\|x^{(1)}-x^{(0)}\|} }{ \ln q } = \dfrac{ \ln \dfrac{10^{-4}(1-\frac{3}{4})}{5} }{ \ln \frac{3}{4} } ≈ 42.43

故根据事前估计，Jacobi约需要迭代43次

Gauss-Seide

{BGS=−(D+L)−1U=[0−2/5−1/50−1/10−11/2001/20−1/8]fGS=(D+L)−1b=[−1252252110]T∥BGS∥∞=1320=0.65<1\begin{cases} B_{GS} = -(D+L)^{-1}U = \left[\begin{array}{c} 0 & −2/5 & −1/5 \\ 0 & −1/10 & −11/20 \\ 0 & 1/20 & −1/8 \end{array}\right] \\ \\ f_{GS} = (D+L)^{-1}b = \left[\begin{array}{c} -\dfrac{12}{5} & \dfrac{22}{5} & \dfrac{21}{10} \end{array}\right]^T \\ \\ \|B_{GS}\|_{\infty} = \dfrac{13}{20} = 0.65 < 1 \end{cases}

x(1)=BGSx(0)+fGS=[−1252252110]Tx^{(1)} = B_{GS}x^{(0)} + f_{GS} = \left[\begin{array}{c} -\dfrac{12}{5} & \dfrac{22}{5} & \dfrac{21}{10} \end{array}\right]^T

∥x(1)−x(0)∥∞=225=4.4\|x^{(1)}-x^{(0)}\|_{\infty} = \dfrac{22}{5} = 4.4

k>ln⁡ε(1−q)∥x(1)−x(0)∥ln⁡q=ln⁡10−4(1−0.65)4.4ln⁡0.65≈27.26k > \dfrac{ \ln \dfrac{ε(1-q)}{\|x^{(1)}-x^{(0)}\|} }{ \ln q } = \dfrac{ \ln \dfrac{10^{-4}(1-0.65)}{4.4} }{ \ln 0.65 } ≈ 27.26

故根据事前估计，Gauss-Seide约需要迭代28次

逐次超松弛迭代（SOR）

①①

记M=1ω(D−ωL)M=\dfrac{1}{ω}(D-ωL)得

B=I−ω(D−ωL)−1A=(D−ωL)−1[(1−ω)D+ωU]B = I - ω(D-ωL)^{-1}A = (D-ωL)^{-1} \big[ (1-ω)D+ωU \big]

②②

记M=1ω(D+ωL)M=\dfrac{1}{ω}(D+ωL)得

B=I−ω(D+ωL)−1A=−(D+ωL)−1[(1−ω)D−ωU]B = I - ω(D+ωL)^{-1}A = -(D+ωL)^{-1} \big[ (1-ω)D-ωU \big]

计算公式

由①①得

Dx(k+1)=Dx(k)+ω(Lx(k+1)+Ux(k)−Dx(k)+b)Dx^{(k+1)} = Dx^{(k)} + ω(Lx^{(k+1)} + Ux^{(k)} - Dx^{(k)} + b)

xi(k+1)=(1−ω)xi(k)+ω⋅bi−∑j=1i−1aijxj(k+1)−∑j=i+1naijxj(k)aiix_i^{(k+1)} = (1-ω)x_i^{(k)} + ω ⋅ \dfrac{ b_i - \sum\limits_{j=1}^{i-1} a_{ij}x_j^{(k+1)} - \sum\limits_{j=i+1}^{n} a_{ij}x_j^{(k)} }{ a_{ii} }

高斯-赛德尔迭代是逐次超松弛迭代的一个特例情况。

ω<1ω<1 ω=1ω=1 ω>1ω>1
低松弛迭代 高斯-赛德尔迭代 超松弛迭代

迭代法的收敛性

定义

误差向量

{e(k+1)=Be(k)e(k)=Bke(0)k∈N\begin{cases} e^{(k+1)} = Be^{(k)} \\ e^{(k)} = B^ke^{(0)} \\ k∈N \end{cases}

收敛性定义

记Ak=(aij(k)∈Rn×n)A_k = (a_{ij}^{(k)}∈R^{n×n})，如果由n2n^2个数列极限存在，且有lim⁡k→∞aij(k)=aij   (i,j=1,2,⋯ ,n)\lim\limits_{k→∞} a_{ij}^{(k)} = a_{ij} \ \ \ (i,j = 1,2,\cdots,n)，则称{Ak}\{A_k\}收敛于AA，记为lim⁡k→∞Ak=A\lim\limits_{k→∞} A_k = A。

收敛性判别

收敛的充要条件

B=(bij)∈Rn×nB=(b_{ij})∈R^{n×n}

lim⁡k→∞Bk=0   ⇔   ρ(B)<1\lim\limits_{k→∞} B^k = 0 { \ \ \ ⇔ \ \ \ } ρ(B)<1

其中ρ(B)ρ(B)为BB的谱半径。

收敛的充分条件

实际计算中，谱半径可能比较难求，可利用ρ(B)≤∥B∥<1ρ(B)≤\|B\|<1（即谱半径≤矩阵的某种算子范数）进行判别

∥B∥=q<1   ⇒   lim⁡k→∞Ak=A\|B\| = q < 1 { \ \ \ ⇒ \ \ \ } \lim\limits_{k→∞} A_k = A

误差估计

qq表示迭代矩阵的某种算子范数

事后估计

若迭代矩阵BB收敛，且允许误差为εε时，当满足

∥x(k)−x(k−1)∥≤1−qqε\|x^{(k)}-x^{(k-1)}\| ≤ \dfrac{1-q}{q}ε

时，有∥x(k)−x(k−1)∥<ε\|x^{(k)}-x^{(k-1)}\| < ε，即可停止迭代。

事前估计

若有

qk1−q∥x(1)−x(0)∥≤ε\dfrac{q^k}{1-q}\|x^{(1)}-x^{(0)}\| ≤ ε

就有∥x(k)−x∗∥<ε\|x^{(k)}-x^{*}\| < ε，即当

k>ln⁡ε(1−q)∥x(1)−x(0)∥ln⁡qk > \dfrac{ \ln \dfrac{ε(1-q)}{\|x^{(1)}-x^{(0)}\|} }{ \ln q }

k∗=⌊k⌋+1k^* = \left\lfloor{ k }\right\rfloor + 1时即可停止迭代。事前估计往往大于实际次数。

收敛速度

由e(k)=Be(k−1)=⋯=Bke(0)e^{(k)} = Be^{(k-1)} = \cdots = B^ke^{(0)}得

∥e(k)∥=∥Bke(0)∥≤∥Bk∥⋅∥e(0)∥\|e^{(k)}\| = \|B^ke^{(0)}\| ≤ \|B^k\| \cdot \|e^{(0)}\|

即

∥Bk∥≥∥e(k)∥∥e(0)∥\|B^k\| ≥ \dfrac{\|e^{(k)}\|}{\|e^{(0)}\|}

平均收敛速度

迭代步数kk与−1kln⁡∥Bk∥-\frac{1}{k}\ln \|B^k\|成反比

Rk(B)=−1kln⁡∥Bk∥R_k(B) = -\dfrac{1}{k}\ln \|B^k\|

迭代矩阵BB的算子范数越小，收敛速度越快。

渐进收敛速度

lim⁡k→∞Rk(B)=−ln⁡(ρ(B))\lim\limits_{k→∞} R_k(B) = -\ln(ρ(B))

R(B)=−ln⁡(ρ(B))R(B) = -\ln(ρ(B))

迭代矩阵BB的谱半径ρ(B)ρ(B)越小，收敛速度越快。

特殊方程收敛性

对角占优阵

设A=(aij)n×nA=(a_{ij})_{n×n}

性质