迭代法解线性方程组

1.迭代法

设有线性方程组
$$\{\begin{array}{llccccccc}{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,& \\ ⋮& & & & & & & & \\ a_{n1}{x}_1& +& a_{n2}{x}_2& +& \cdots & +a_{nn}{x}_n& =& b_n,& \end{array}$$现将该线性方程组改写为
$$\{\begin{array}{llcccc}{x}_1=\frac{1}{a_{11}}(-a_{12}{x}_2& -\cdots & & -a_{1n}{x}_n& +& b_1),\\ x_2=\frac{1}{a_{22}}(-a_{21}{x}_1& -a_{23}{x}_3& -\cdots & -a_{2n}{x}_n& +& b_2),\\ ⋮& & & & & \\ x_n=\frac{1}{a_{nn}}(-a_{n1}{x}_1& -a_{n2}{x}_2& -\cdots & -a_{nn-1}{x}_{n-1}& +& b_n);(1)\end{array}$$或者写为$x=B_{0}x+f,$其中
$$B_0=\left[\begin{array}{ccccc}0& -\frac{a_{12}}{a_{11}}& -\frac{a_{13}}{a_{11}}& \dots & -\frac{a_{1n}}{a_{11}}\\ -\frac{a_{21}}{a_{22}}& 0& -\frac{a_{23}}{a_{22}}& \dots & -\frac{a_{2n}}{a_{22}}\\ -\frac{a_{31}}{a_{33}}& -\frac{a_{32}}{a_{33}}& 0& \dots & -\frac{a_{3n}}{a_{33}}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ -\frac{a_{n1}}{a_{nn}}& -\frac{a_{n2}}{a_{nn}}& \dots & -\frac{a_{nn-1}}{a_{nn}}& 0\end{array}\right],f=\left[\begin{array}{c}\frac{b_1}{a_{11}}\\ \frac{b_2}{a_{22}}\\ \frac{b_3}{a_{33}}\\ ⋮\\ \frac{b_n}{a_{nn}}\end{array}\right]$$
任取初始值$x^{(0)}=(0,0,\dots ,0{)}^T$,将其带入(1)式右边，得到新的值$x^{(1)}=(x_1^{(1)},x_2^{(1)},\dots ,x_n^{(1)}{)}^T$
再将$x^{(1)}$带入右边，反复进行从而可以得到一个向量序列
$$x^{(k+1)}=B{x}^{(k)}+f,k=0,1,2,\dots ,(2)$$其中k表示迭代次数.
迭代法及其收敛性定义
(1)对于给定的线性方程组$x=Bx+f,$用公式(2)逐步带入求近似解的方法称为迭代法(或一阶定常迭代法,这里$B与k$无关).
(2)如果$\underset{k\to \infty }{\lim }{x}^{(k)}存在(记为x^{\ast }),$称次迭代法收敛，显然$x^{\ast }$就是此方程的解，否则称此迭代法发散.

2.雅可比迭代法(Jacobi)

将线性方程组$Ax=b$中的系数矩阵$A=(a_{ij})\in {\mathbb{R}}^{n\times n}$分成三部分
$$A=\left[\begin{array}{cccc}{a}_{11}& & & \\ & a_{22}& & \\ & & \ddots & \\ & & & a_{nn}\end{array}\right]-\left[\begin{array}{ccccc}0& & & & \\ -a_{21}& 0& & & \\ ⋮& ⋮& \ddots & & \\ -a_{n-1,1}& -a_{n-1,2}& \dots & 0& \\ -a_{n1}& -a_{n2}& \dots & -a_{n,n-1}& 0\end{array}\right]-\left[\begin{array}{ccccc}0& -a_{12}& \dots & -a_{1,n-1}& -a_{1n}\\ & 0& \dots & -a_{2,n-1}& -a_{2n}\\ & & \ddots & ⋮& ⋮\\ & & & 0& -a_{n-1,n}\\ & & & & 0\end{array}\right]\equiv D-L-U.$$
设$a_{ii}\ne 0(i=1,2,\dots ,n),$选取$M$为$A$的对角元素部分，即选取$M=D(对角矩阵),A=D-N,$得到解$Ax=b$的雅可比(Jacobi)迭代法
$$\{\begin{array}{l}{x}^{(0)},\text{初始向量,}\\ x^{(k+1)}=B{x}^{(k)}+f,k=0,1,2,\dots ,\end{array}$$
其中$B=I-D^{-1}A=D^{-1}(L+U)\equiv J,f=D^{-1}b.$称$J为Ax=b$的雅可比迭代法的迭代矩阵.

3.高斯-赛德尔迭代法(Gauss-Seidel)

选取分裂矩阵$M为A$的下三角部分，即选取$M=D-L(下三角矩阵),A=M-N,于是得到解Ax=b的高斯-赛德尔迭代法.$
$$\{\begin{array}{l}{x}^{(0)},\text{初始向量,}\\ x^{(k+1)}=B{x}^{(k)}+f,k=0,1,2,\dots ,\end{array}$$
其中$B=I-(D-L{)}^{-1}A=(D-L{)}^{-1}U\equiv G,f=(D-L{)}^{-1}b.$
称$G=(D-L{)}^{-1}U为解Ax=b$的高斯-赛德尔迭代法的迭代矩阵.

4.逐次超松弛迭代法(SOR)

选取分裂矩阵$M$为带参数的下三角矩阵$$M=\frac{1}{\omega }(D-\omega L),$$
其中$\omega >0$为可选择的松弛因子.于是可以构造一个迭代矩阵
$$L_{\omega }\equiv I-\omega (D-\omega L{)}^{-1}A=(D-\omega L{)}^{-1}((1-\omega )D+\omega U).$$
从而得到解$Ax=b$的逐次超松弛迭代法(successive over relaxation method,简称SOR方法).
解$Ax=b$的SOR方法为
$$\{\begin{array}{l}{x}^{(0)},\text{初始向量,}\\ x^{(k+1)}=L_{\omega }{x}^{(k)}+f,k=0,1,\dots ,\end{array}$$
其中$L_{\omega }=(D-\omega L{)}^{-1}((1-\omega )D+\omega U),f=\omega (D-\omega L{)}^{-1}b.$

5.Python实现三种迭代法

import numpy as np
from scipy.linalg import eigvals, inv

def vector_norm(vector: np.ndarray, p=None):
    """
    计算向量的p-范数
    :param vector: 实向量或者复向量
    :param p: 指定类型的范数，默认是oo范数
    :return: 指定向量范数和p值
    """
    if p is None:
        return abs(vector).max(), p
    elif p >= 1:
        return np.power(np.sum(np.power(abs(vector), p)), 1 / p), p
    else:
        raise Exception("error,p must be an integer , greater than  or equal to 1")


def jacobi(coefficient_matrix: np.ndarray,
           right_hand_side_vector: np.ndarray,
           initial_vector: np.ndarray, epsilon=1e-4):
    d = np.diag(np.diag(coefficient_matrix))
    l = -np.tril(coefficient_matrix, k=-1)
    u = -np.triu(coefficient_matrix, k=1)
    jacobi_matrix = inv(d) @ (l + u)
    if spectral_radius(jacobi_matrix) < 1:
        f = inv(d) @ right_hand_side_vector
        x_1 = initial_vector
        x_2 = jacobi_matrix @ x_1 + f
        iteration_count = 1
        while 1:
            norm, _ = vector_norm(x_2 - x_1)
            if norm < epsilon:
                return x_2, iteration_count
            x_1 = x_2
            x_2 = jacobi_matrix @ x_1 + f
            iteration_count += 1
    raise Exception("jacobi iteration is not convergent")


def gauss_seidel(coefficient_matrix: np.ndarray,
                 right_hand_side_vector: np.ndarray,
                 initial_vector: np.ndarray, epsilon=1e-4):
    d = np.diag(np.diag(coefficient_matrix))
    l = -np.tril(coefficient_matrix, k=-1)
    u = -np.triu(coefficient_matrix, k=1)
    gauss_seidel_matrix = inv(d - l) @ u
    if spectral_radius(gauss_seidel_matrix) < 1:
        f = inv(d - l) @ right_hand_side_vector
        x_1 = initial_vector
        x_2 = gauss_seidel_matrix @ x_1 + f
        iteration_count = 1
        while 1:
            norm, _ = vector_norm(x_2 - x_1)
            if norm < epsilon:
                return x_2, iteration_count
            x_1 = x_2
            x_2 = gauss_seidel_matrix @ x_1 + f
            iteration_count += 1
    raise Exception("jacobi iteration is not convergent")


def successive_over_relaxation(coefficient_matrix: np.ndarray,
                               right_hand_side_vector: np.ndarray,
                               initial_vector: np.ndarray,
                               true_solution: np.ndarray, omega=1.0, epsilon=5e-6):
    d = np.diag(np.diag(coefficient_matrix))
    l = -np.tril(coefficient_matrix, k=-1)
    u = -np.triu(coefficient_matrix, k=1)
    l_omega = inv(d - omega * l) @ ((1 - omega) * d + omega * u)
    if spectral_radius(l_omega) < 1:
        f = omega * inv(d - omega * l) @ right_hand_side_vector
        x = l_omega @ initial_vector + f
        iteration_count = 1
        while 1:
            norm, _ = vector_norm(true_solution - x)
            if norm < epsilon:
                return x, iteration_count
            x = l_omega @ x + f
            iteration_count += 1
    raise Exception("jacobi iteration is not convergent")


def hilbert_matrix(order: int):
    hilbert = np.zeros((order, order))
    for i in range(order):
        for j in range(order):
            hilbert[i, j] = 1 / (i + j + 1)
    return hilbert


def spectral_radius(square_matrix: np.ndarray):
    if square_matrix.shape[0] == square_matrix.shape[1]:
        return abs(eigvals(square_matrix)).max()
    raise Exception("\n{} is not a square matrix".format(square_matrix))