数值计算方法 Chapter6. 解线性方程组的迭代法

codename_cys
230 阅读5分钟

本文已参与「新人创作礼」活动，一起开启掘金创作之路。

0. 问题描述

这一章节要解的问题和上一章是一样的，依然还是nn元线性方程组的求解问题。

{a11x1+a12x2+...+a1nxn=y1a21x1+a22x2+...+a2nxn=y2...an1x1+an2x2+...+annxn=yn\left\{ \begin{aligned} a_{11}x_1 + a_{12}x_2 + ... + a_{1n}x_n &= y_1 \\ a_{21}x_1 + a_{22}x_2 + ... + a_{2n}x_n &= y_2 \\ ... \\ a_{n1}x_1 + a_{n2}x_2 + ... + a_{nn}x_n &= y_n \end{aligned} \right.

但是，上一章主要是通过矩阵的线性变换转换成可以快速求解的三角阵或者对角阵的方式进行求解，其计算结果是精确的结果。

而本章则是的思路则是将问题Ax=yAx=y转换成x=A′xx = A'x的迭代形式，从而，我们就可以给出迭代数组xk+1=A′xkx_{k+1} = A'x_{k}。

此时，如果xkx_k满足收敛条件，那么xkx_{k}就会收敛到x=A′xx = A'x的一组解当中，上述问题同样可以得到解答。

1. Jacobi迭代

1. Jacobi迭代方法

对原始方程进行线性变换，我们显然有：

{x1=−a12a11x2−...−a1na11xn+y1a11x2=−a21a22x2−...−a2na22xn+y2a22...xn=−an1annx1−an2annx2−...+ynann\left\{ \begin{aligned} x_{1} = & &-\frac{a_{12}}{a_{11}}x_2 &-... &-\frac{a_{1n}}{a_{11}}x_n &+ \frac{y_1}{a_{11}} \\ x_{2} = &-\frac{a_{21}}{a_{22}}x_2 & &- ... &-\frac{a_{2n}}{a_{22}}x_n &+ \frac{y_2}{a_{22}} \\ ... \\ x_{n} = &-\frac{a_{n1}}{a_{nn}}x_1 &-\frac{a_{n2}}{a_{nn}}x_2 &- ... & &+ \frac{y_n}{a_{nn}} \end{aligned} \right.

我们将其写成矩阵形式即为：

(x1x2...xn)=(0−a12a11...−a1na11−a21a220...−a2na22............−an1ann−an2ann...0)(x1x2...xn)+(y1a11y2a22...ynann)\begin{pmatrix} x_1 \\ x_2 \\ ... \\ x_n \end{pmatrix} = \begin{pmatrix} 0 & -\frac{a_{12}}{a_{11}} & ... & -\frac{a_{1n}}{a_{11}} \\ -\frac{a_{21}}{a_{22}} & 0 & ... & -\frac{a_{2n}}{a_{22}} \\ ... & ... & ... & ... \\ -\frac{a_{n1}}{a_{nn}} & -\frac{a_{n2}}{a_{nn}} & ... & 0 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ ... \\ x_n \end{pmatrix} + \begin{pmatrix} \frac{y_1}{a_{11}} \\ \frac{y_2}{a_{22}} \\ ... \\ \frac{y_n}{a_{nn}} \end{pmatrix}

我们将其简化为符号表达即为：

x=Bx+gx = Bx + g

基于此，我们就可以给出Jacobi迭代公式如下：

xk+1=Bxk+gx_{k+1} = Bx_{k} + g

2. Jacobi迭代矩阵

我们给出更加严格的数据表达如下，令：

D=(a11a22...ann)D = \begin{pmatrix} a_{11} \\ & a_{22} \\ & & ... \\ & & & a_{nn} \end{pmatrix}

则有：

Ax=(D+A−D)x=yAx = (D + A - D)x = y

进而可以导出：

Dx=(D−A)x+yDx = (D-A)x + y

两边左乘D−1D^{-1}即有：

x=D−1(D−A)x+D−1yx = D^{-1}(D-A)x + D^{-1}y

记B=D−1(D−A)B=D^{-1}(D-A)，g=D−1yg=D^{-1}y，即可得到上一小节中一模一样的内容。

迭代矩阵BB即为Jacobi迭代矩阵。

3. Jacobi迭代收敛条件

Jacobi迭代收敛的充要条件为：

迭代矩阵BB的谱半径ρ(B)<1\rho(B) < 1。

亦即：

  • 对于迭代矩阵BB的所有本征值λi\lambda_i，均有∣λi∣<1|\lambda_i| < 1，或者等价的max(∣λi∣)<1max(|\lambda_i|) < 1。

不过，在一些特定的情况下，上述充要条件可以有适当的放宽。

具体而言，有如下定理：

定理 6.1
若方程组Ax=yAx=y的系数矩阵AA，满足下列条件之一，则其Jacobi迭代收敛：
(1) A为行对角优阵，即∣aii∣>∑j≠i∣aij∣|a_{ii}| > \sum_{j \neq i} |a_{ij}|，i=1,2,...,ni=1,2,...,n；
(2) A为列对角优阵，即∣ajj∣>∑i≠j∣aij∣|a_{jj}| > \sum_{i \neq j} |a_{ij}|，j=1,2,...,nj=1,2,...,n；

4. python伪代码实现

最后，我们给出Jacobi迭代的python伪代码如下：

def jacobi_iter(A, y, epsilon=1e-6):
    n = len(A)
    B = [[0 for _ in range(n)] for _ in range(n)]
    g = [0 for _ in range(n)]
    for i in range(n):
        for j in range(n):
            if j == i:
                continue
            B[i][j] = -A[i][j] / A[i][i]
        g[i] = y[i] / A[i][i]
    
    x = [0 for _ in range(n)]
    for _ in range(10**6):
        z = [0 for _ in range(n)]
        for i in range(n):
            z[i] += g[i]
            for j in range(n):
                z[i] += B[i][j] * x[j]
        if max(abs(z[i] - x[i]) for i in range(n)) < epsilon:
            break
        x = z
    return z

2. Gauss-Seidel迭代

1. Gauss-Seidel迭代方法

Gauss-Seidel迭代方程和上述Jacobi迭代事实上是非常相似的，唯一的区别在于说Jacobi迭代是以x(k)x^{(k)}为整体每次一起进行迭代更新的，而Guass-Seidel迭代则是在计算每一个xi(k)x^{(k)}_{i}的时候就是用当前已经迭代计算完成的所有的xj(k)x^{(k)}_{j}的结果。

即是说，以每一个元素为单位不断进行迭代更新。

用公式表达即为：

{x1(k+1)=−a12a11x2(k)−...−a1na11xn(k)+y1a11x2(k+1)=−a21a22x1(k+1)−...−a2na22xn(k)+y2a22...xn(k+1)=−an1annx1(k+1)−an2annx2(k+1)−...+ynann\left\{ \begin{aligned} x^{(k+1)}_{1} = & &-\frac{a_{12}}{a_{11}}x^{(k)}_2 &-... &-\frac{a_{1n}}{a_{11}}x^{(k)}_n &+ \frac{y_1}{a_{11}} \\ x^{(k+1)}_{2} = &-\frac{a_{21}}{a_{22}}x^{(k+1)}_1 & &- ... &-\frac{a_{2n}}{a_{22}}x^{(k)}_n &+ \frac{y_2}{a_{22}} \\ ... \\ x^{(k+1)}_{n} = &-\frac{a_{n1}}{a_{nn}}x^{(k+1)}_1 &-\frac{a_{n2}}{a_{nn}}x^{(k+1)}_2 &- ... & &+ \frac{y_n}{a_{nn}} \end{aligned} \right.

2. Gauss-Seidel迭代矩阵

同样的，我们给出更加严格的数学推导如下：

A=D+L+U=(a11a22...ann)+(0a210...an1an2...0)(0a12...a1n0...a2n...0)\begin{aligned} A &= D + L + U \\ &= \begin{pmatrix} a_{11} \\ & a_{22} \\ & & ... \\ & & & a_{nn} \end{pmatrix} + \begin{pmatrix} 0 \\ a_{21} & 0 \\ ... \\ a_{n1} & a_{n2} & ... & 0 \end{pmatrix} \begin{pmatrix} 0 & a_{12} & ... & a_{1n} \\ & 0 & ... & a_{2n} \\ & & ... \\ & & & 0 \end{pmatrix} \end{aligned}

从而，我们有：

Ax=(D+L+U)x=yAx = (D+L+U)x = y

亦即：

(D+L)x=−Ux+y(D+L)x = -Ux + y

两边同乘以(D+l)−1(D+l)^{-1}即有：

x=−(D+L)−1Ux+(D+L)−1yx = -(D+L)^{-1}Ux + (D+L)^{-1}y

令S=−(D+L)−1US = -(D+L)^{-1}U，f=(D+L)−1yf = (D+L)^{-1}y，我们即可写出Gauss-Seidel迭代公式：

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

称迭代矩阵SS即为Gauss-Seidel迭代矩阵。

当然，如上一小节所述，实际在算法实现当中并没有必要计算出SS和ff，只需要根据前面JacobiJacobi迭代矩阵进行实时参数更新即可。

3. Gauss-Seidel迭代收敛条件

同样的，我们给出书中关于Gauss-Seidel迭代的收敛条件如下：

定理6.2
若方程组系数矩阵为行或列对角优时，则Gauss-Seidel迭代收敛。
定理6.3
若方程组系数矩阵AA为对称正定阵，则Gauss-Seidel迭代收敛。

4. 伪代码实现

同样的，我们用python给出伪代码如下：

def gauss_seidel_iter(A, y, epsilon=1e-6):
    n = len(A)
    B = [[0 for _ in range(n)] for _ in range(n)]
    g = [0 for _ in range(n)]
    for i in range(n):
        for j in range(n):
            if j == i:
                continue
            B[i][j] = -A[i][j] / A[i][i]
        g[i] = y[i] / A[i][i]
    
    x = [0 for _ in range(n)]
    cnt = 0
    for _ in range(10**6):
        for i in range(n):
            z = g[i]
            for j in range(n):
                z += B[i][j] * x[j]
            if abs(z-x[i]) < epsilon:
                cnt += 1
            else:
                cnt = 0
            x[i] = z

            if cnt >= n:
                break
        if cnt >= n:
            break
    return x

3. 松弛迭代

1. 松弛迭代方法

松弛迭代的原型依然还是之前的Jacobi迭代，不过，和Gauss-Seidel迭代的实时参数更新不同，松弛迭代在这里是对Jacobi迭代式的批次更新以及Gauss-Seidel迭代式的实时更新取了一个折中，通过一个超参ww来进行参数更新速度的控制。

具体而言，即为：

{x1(k+1)=(1−w)x1(k)+w(b12x2(k)+b13x3(k)+...+b1nxn(k)+g1)x2(k+1)=(1−w)x2(k)+w(b21x1(k+1)+b23x3(k)+...+b2nxn(k)+g2)...xn(k+1)=(1−w)xn(k)+w(bn1x1(k+1)+bn2x2(k+1)+...+bn,n−1xn−1(k+1)+gn)\left\{ \begin{aligned} x^{(k+1)}_1 &= (1-w)x^{(k)}_1 + w(b_{12}x^{(k)}_2 + b_{13}x^{(k)}_3 + ... + b_{1n}x^{(k)}_n + g_1) \\ x^{(k+1)}_2 &= (1-w)x^{(k)}_2 + w(b_{21}x^{(k+1)}_1 + b_{23}x^{(k)}_3 + ... + b_{2n}x^{(k)}_n + g_2) \\ ... \\ x^{(k+1)}_n &= (1-w)x^{(k)}_n + w(b_{n1}x^{(k+1)}_1 + b_{n2}x^{(k+1)}_2 + ... + b_{n,n-1}x^{(k+1)}_{n-1} + g_n) \end{aligned} \right.

2. 松弛迭代矩阵

同样的，我们给出松弛迭代的数学表达如下：

x(k+1)=(1−w)x(k)+w(L~x(k+1)+U~x(k)+g)x^{(k+1)} = (1-w)x^{(k)} + w(\tilde{L}x^{(k+1)} + \tilde{U}x^{(k)} + g)

仿照Gauss-Seidel迭代，我们同样可以给出其迭代矩阵格式如下：

x(k+1)=Swx(k)+fx^{(k+1)} = S_w x^{(k)} + f

其中，

{Sw=(I+wD−1L)−1[(1−w)I−wD−1U]f=w(I+wD−1L)−1g\left\{ \begin{aligned} S_w &= (I + wD^{-1}L)^{-1}[(1-w)I - wD^{-1}U] \\ f &= w(I + wD^{-1}L)^{-1}g \end{aligned} \right.

3. 松弛迭代的收敛条件

而松弛迭代的收敛判断存在如下两个定理：

定理 6.4
松弛迭代收敛的必要条件为0<w<20 < w < 2；
定理 6.5
若AA为对称正定矩阵，则当0<w<20 < w < 2时，松弛迭代恒收敛；

特别的，当0<w<10 < w < 1时，称其为亚松弛迭代；当1<w<21 < w < 2时，称其为超松弛迭代；当w=1w=1时，迭代退化为Gauss-Seidel迭代。

4. 伪代码实现

最后，我们同样给出松弛迭代的伪代码实现如下：

def loose_iter(A, y, w=0.5, epsilon=1e-6):
    n = len(A)
    B = [[0 for _ in range(n)] for _ in range(n)]
    g = [0 for _ in range(n)]
    for i in range(n):
        for j in range(n):
            if j == i:
                continue
            B[i][j] = -A[i][j] / A[i][i]
        g[i] = y[i] / A[i][i]
    
    x = [0 for _ in range(n)]
    cnt = 0
    for _ in range(10**6):
        for i in range(n):
            z = (1-w) * x[i] + w * g[i]
            for j in range(n):
                z += w * B[i][j] * x[j]
            if abs(z-x[i]) < epsilon:
                cnt += 1
            else:
                cnt = 0
            x[i] = z

            if cnt >= n:
                break
        if cnt >= n:
            break
    return x

4. 逆矩阵计算

最后，我们来考察一下逆矩阵计算。

逆矩阵的计算原则上来说其实算是上述解线性方程组的一个特殊应用，事实上解nn个单元向量然后将其解拼接一下就能得到我们的逆矩阵了。

我们这里就不进行详述了，只给出python伪代码实现如下：

def cal_inverse_matrix(A):
    n = len(A)
    res = [[0 for _ in range(n)] for _ in range(n)]
    for j in range(n):
        y = [0 for _ in range(n)]
        y[j] = 1
        x = gauss_seidel_iter(A, y)
        for i in range(n):
            res[i][j] = x[i]
    return res
avatar
文章
阅读
粉丝