对角占优线性方程组的迭代法

线性方程组的求解是多数数值问题最终的一个环节，很多问题最终都变为了一个求解线性方程组的问题。求解方程组最广为人知的方法是高斯消元法，这是一种直接方法，通过构造某种形式，可以直接得到最终的结果，其矩阵形式就对应着LU分解法。现代的计算机语言对于科学计算的解决，其已经包含了直接法，并且实现的效率极高，而迭代法计算效率低下，一般计算中往往不会采用它。但是，迭代法求解问题的思路却是十分有价值的，迭代思想也不仅仅局限于方程组的求解。

迭代的直观感受

以一个二元一次方程组为例，

$$ a_{11}x_{1} + a_{12}x_{2} = b_{1}\\ a_{21}x_{1} + a_{22}x_{2} = b_{2} $$

从方程组的第一式得到：

$$ x_{1} = \frac{b_{1} - a_{12}x_{2}}{a_{11}} $$

第二式得到：

$$ x_{2} = \frac{b_{2} - a_{21}x_{1}}{a_{21}} $$

以上构造的关系式，再附加上一组初始值，通过反复的计算就可以得到方程组的近似解。上述的迭代公式的构造称为不动点迭代，该院里在方程求解中有过应用，这里该方法也称为雅克比法，Jacobian法。

Jacobian法

前述的实例可以看出，jacobian法本质上就是对每一个方程进行变换，依次得到一个$x_i$的表达式，然后依据初始值就可以玩车呢个迭代计算，且这一方法中，每一个式均是对应$b$减去其他项的和，然后除以$x_i$的系数。这样可以方便地得到其矩阵形式。

方程组为$Ax = b$，对系数矩阵$A$进行分解，将其分解为对角线一下元素组成的矩阵$L$，只有对角线元素的矩阵$D$，以及对角线元素以上的元素组成的矩阵$U$，这样矩阵$A = L + D + U$，由此得到：

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

得到一个迭代公式：

$$ X_{k + 1} = D^{-1} b - D^{-1}(L + U) X_{k} \quad (k = 0,1,2,3,....) $$

这样第一次的给定的初始值记做$X_{0}$，利用它可以得到一组近似值为$X_{1}$，反复迭代直到精度达到事先设定的程度。这一方法可以很方便地利用代码实现。

import numpy as np
# 定义一个函数来完成jacobian法
# 输入参数包括方程组本身、初始值、精度
# 这里为了避免陷入死循环，还需要一个控制最大迭代次数

def Jacobian_Method(A, b, x0, eps=0.0001, max_times=400):

    L = np.tril(A, -1)
    U = np.triu(A, 1)
    D = np.diag(np.diag(A, 0))
    number= 0  # 记录迭代次数
    while np.linalg.norm(np.dot(A, x0) - b, np.inf) > eps:
        F = np.linalg.solve(D, b)
        G = np.linalg.solve(D, L + U)
        x0 = F - np.dot(G, x0)
        number += 1  #此时，迭代完成一次，则次数加一
        if number > max_times:
            print('Warning! The solution is not stable!')
            return
    return x0, number

这是雅克比法的数值实现代码。该方法的收敛速度并不快，并且并非对任意的线性方程组都能够进行求解，迭代法的适用范围主要是针对对角占优矩阵的。

Gauss-Seidel迭代

从Jcobian方法中可以看出的是，每一次的计算中，第一个未知数必定会被先得到确定，那么接下来第二个未知数的计算中人然会用第一个未知数，只不过对于Jcobian来说，它的计算依据是采用的上一步的值，而Gauss-Seidel法采用了更新的值用于后面的未知数的计算，这样可以大大提高收敛的速度。其矩阵形式为：

$$ (L + D + U)x = b\\ （L + D）x = b - Ux\\ x = (L + D)^{-1}B - (L + D)^{-1} U x $$

迭代公式为：

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

def Gauss_Seidel_Method(A, b, x0, eps=0.0001, max_times=400):
    L = np.tril(A, -1)
    U = np.triu(A, 1)
    D = np.diag(np.diag(A, 0))
    number= 0  # 记录迭代次数
    while np.linalg.norm(np.dot(A, x0) - b, np.inf) > eps:
        F = np.linalg.solve(L + D, b)
        G = np.linalg.solve(L + D, U)
        x0 = F - np.dot(G, x0)
        number += 1
        if number > max_times:
            print('Warning: The result is not stable!')
            return
    return x0, number

Gauss-Seidel迭代法的收敛速度也明显快于Jacobian法，所以这一方法使用的更为广泛。但是，有没有收敛速度更快的方法呢？这就是连续松弛迭代法，它本质是对Gauss-Seidel法的一个改进，通过松弛参数$\omega$。

连续松弛迭代法SOR

在gauss-seidel法的基础上，通过对$Ax=b$两边乘以一个松弛参数$\omega$，再进行迭代公式的构建。

$$ \omega A x = \omega b\\ \omega (L + D + U) X = \omega b\\ \omega Lx + \omega Dx + \omega Ux + Dx = \omega b + Dx\\ (\omega L + D)x = \omega b + Dx - \omega D x - \omega U x\\ x = (\omega L + D)^{-1}[(1 - \omega)D - \omega U]x + \omega (\omega L + D)^{-1} b $$

得到迭代公式

$$ X_{k+1} = (\omega L + D)^{-1}[(1 - \omega)D - \omega U]X_{k} + \omega (\omega L + D)^{-1} b $$

这一方法的收敛速度与$\omega$有关，当大于1时为过松弛迭代，小于1时为欠松弛迭代，等于1是就完全退化为gauss-seidel法。而且对于松弛参数的确定并没有一个确定的方法，需要通过不断的试算采纳呢个估计到一个最优的值。

def SOR_Method(A, b, x0, w, eps=0.0001, max_times=400):
    L = np.tril(A, -1)
    U = np.triu(A, 1)
    D = np.diag(np.diag(A, 0))
    number= 0  # 记录迭代次数
    while np.linalg.norm(np.dot(A, x0) - b, np.inf) > eps:
        F = w * np.linalg.solve((w * L + D), b)
        G = np.linalg.solve(w * L + D, (1 - w) * D - w * U)
        x0 = F + np.dot(G, x0)
        number += 1
        if number > max_times:
            print('Warning: The result is not stable!')
            return
    return x0, number

收敛条件与效率对比

上述的迭代法是否一定收敛？这个会收到系数矩阵$A$的影响，如果系数矩阵一个严格的对角占优矩阵，那么上述的三个迭代法一定收敛，但是收敛并不意味着这个矩阵一定是严格的对角占优矩阵，严格对角占优只是收敛的充分条件。对角占优表示每一行的主元比该行其他所有元素之和的绝对值还要大，这就是对角占优。

在线性方程组的求解中，迭代法现在一般不会被使用，主要原因是迭代法的计算效率差，没有直接法快速，计算机硬件的普及，使得较大规模的直接发求解很容易实现，并且数值计算语言中也已经设置好了算法，matlab中的"\"就是直接法，这是一个高级的LU分解法的实现，并且利用底层的c++进行了预编译，速度十分快，而自己定义的迭代法在没有编译的情况下，速度慢大约近100倍，在python中，numpy库实现了线性代数的基本运算，这里的solve方法相当于matlab的"\"运算。

A = np.array([[3,1,-1],[2,4,1],[-1,2,5]])
b= np.array([4,1,1])
x0 = [0,0,0]
w = 1.1
a1 = Jacobian_Method(A, b, x0)
a2 = Gauss_Method(A, b, x0)
a3 = SOR_Method(A, b, x0, w)
print(a1)
print(a2)
print(a3)

结果如下：

(array([ 1.99995   , -0.99994562,  0.99995301]), 26)
(array([ 1.99996414, -0.99997146,  0.99998141]), 13)
(array([ 1.9999659 , -0.99997676,  0.9999868 ]), 9)

雅克比法迭代了26次才达到精度，高斯赛德尔法用了13次，SOR中，当松弛参数为1.1时，连续过松弛方法的收敛用了9次。

本站文章如非特别说明，均为原创。未经本人同意，请勿转载！转载请务必注明出处！