Python实现Gauss-Seidel迭代法-CFANZ编程社区

在雅可比（Jacobi）迭代法的基础上，通过修改其中逻辑使其成为Gauss-seidel迭代法。通过输入系数矩阵mx，值矩阵mr，最大迭代次数n，目标误差e即可得到答案。

代码如下：

import numpy as np


def Astringency(mx):  # 判断系数矩阵的收敛性
    L, D, U = [], [], []  # 初始化L，D，U矩阵
    for i in range(len(mx)):
        L.append([]), D.append([]), U.append([])
        for j in range(len(mx)):
            if i > j:
                L[i].append(mx[i][j]), D[i].append(0), U[i].append(0)
            if i == j:
                L[i].append(0), D[i].append(mx[i][j]), U[i].append(0)
            if i < j:
                L[i].append(0), D[i].append(0), U[i].append(mx[i][j])
    # print(L)
    # print(D)
    # print(U)
    ld = L  # 计算D+L矩阵
    for i in range(len(mx)):
        for k in range(len(mx)):
            ld[i][k] = L[i][k] + D[i][k]
    # print(ld)
    # print(-np.linalg.inv(ld))
    # print(D)
    G = np.dot(-np.linalg.inv(ld), U)  # 得到G矩阵
    # print(G)
    e, v = np.linalg.eig(G)
    # print(e)
    for i in range(len(e)):
        count = 0  # 计算不小于一的特征值的数量
        if abs(e[i]) < 1:
            continue
        else:
            count = count + 1
            print(abs(e[i]))
    # print(count)
    if count == 0:
        return True
    else:
        print("迭代法不收敛")
        return False


def Gauss_seidel(mx, mr, n=100, e=0.0001):  # mx为系数矩阵，mr为值矩阵,n为默认迭代次数，e为默认误差返回
    if len(mx) == len(mr):  # 若mx与mr长度相等则开始迭代，否则方程无解
        if Astringency(mx) == 1:  # 判断系数矩阵mx是否收敛
            x = []
            for i in range(len(mr)):
                x.append([0])  # 得到长度与mr相等的初值，并且把初值设定为零
            count = 0  # 迭代次数计数
            while count < n:  # 没有达到迭代次数时继续迭代
                tempx = np.copy(x)  # tempx用于存储迭代前的x值,但又不指向同一对象，用于后面的误差计算
                for i in range(len(x)):
                    ri = mr[i][0]
                    for k in range(len(mx[i])):
                        if k != i:
                            ri = ri - mx[i][k] * x[k][0]
                    ri = ri / mx[i][i]
                    x[i][0] = ri  # 每次计算存储单个x值
                print("第{}次迭代的值为：{}".format(count + 1, x))
                ee = []  # 存储每两次迭代结果之间的误差
                for i in range(len(x)):
                    # print(x[i][0])
                    # print(tempx[i][0])
                    ee.append(abs(x[i][0] - tempx[i][0]))
                em = max(ee)  # 取最大误差值
                print("第{}、{}次迭代间误差值为:{}".format(count, count + 1, em))
                if em < e:
                    return x  # 当两次迭代的x的最大误差满足要求时，直接返回计算结果
                count += 1
            return False  # 当运行到最大迭代次数时精度仍不满足要求则返回错误

        else:
            print("使用迭代法不收敛")
    else:
        print("此方程无解")


mx = [[5, 2, 1], [-1, 4, 2], [2, -5, 10]]
# # print(len(mx))

mr = [[-12], [10], [1]]
# print(len(mr))
print(Gauss_seidel(mx, mr, 100, 0.0001))

注意点：存储前一次迭代解时，若采用tempx=x，tempx与x将指向同一对象，会使迭代无法往下进行。解决方法：采用tempx=np.copy(x)创建唯一副本。

实现结果：

第1次迭代的值为：[[-2.4], [1.9], [1.53]]
第0、1次迭代间误差值为:2.4
第2次迭代的值为：[[-3.466], [0.8684999999999999], [1.22745]]
第1、2次迭代间误差值为:1.0660000000000003
第3次迭代的值为：[[-2.99289], [1.1380525000000001], [1.2676042500000002]]
第2、3次迭代间误差值为:0.47311000000000014
第4次迭代的值为：[[-3.1087418500000004], [1.0890124124999998], [1.26625457625]]
第3、4次迭代间误差值为:0.11585185000000031
第5次迭代的值为：[[-3.08885588025], [1.0946587418125002], [1.2651005469562502]]
第4、5次迭代间误差值为:0.019885969750000232
第6次迭代的值为：[[-3.09088360611625], [1.0947288249928122], [1.265541133719656]]
第5、6次迭代间误差值为:0.0020277258662497744
第7次迭代的值为：[[-3.090999756741056], [1.094479493954908], [1.2654396983256653]]
第6、7次迭代间误差值为:0.00024933103790414357
第8次迭代的值为：[[-3.090879737247096], [1.0945602165253931], [1.2654560557121157]]
第7、8次迭代间误差值为:0.00012001949396012179
第9次迭代的值为：[[-3.0909152977525802], [1.094543147705797], [1.2654546334034145]]
第8、9次迭代间误差值为:3.55605054842556e-05
[[-3.0909152977525802], [1.094543147705797], [1.2654546334034145]]

Jacobi迭代法与Gauss-Seidel迭代法比较：
仅此例子比较，Jacobi迭代法进行二十次迭代后达到目标精度，而Gauss-Seidel迭代法仅经过十次迭代后就到达了目标精度。从该例子比较Gauss-Seidel迭代法迭代速度更加快，可能的原因是Gauss-Seidel迭代法的计算数值处于一直更新迭代状态，因此迭代速度较快。