基于python通过迭代公式的矩阵表示实现高斯-赛德尔迭代方法

3.2 线性方程组的求解 (100分)

分别使用迭代法（Jacobi与Gausss-Seidel迭代法）和直接求解法（消元法与三角分解法）求解线性方程组：

要求：使误差不超过。

①请用Jacobi迭代法、Gausss-Seidel迭代法实现，并比较不同迭代方法的计算工作量和计算结果。

②选择消元法与三角分解法中的一种方法实现。

③比较迭代法和直接求解法的异同。

import numpy as np
a=np.mat([[7.2,2.3,-4.4,0.5],
          [1.3,6.3,-3.5,2.8],
          [5.6,0.9,8.1,-1.3],
          [1.5,0.4,3.7,5.9]])
L=np.mat([[0,0,0,0],
          [1.3,0,0,0],
          [5.6,0.9,0,0],
          [1.5,0.4,3.7,0]])
D=np.mat([[7.2,0,0,0],
          [0,6.3,0,0],
          [0,0,8.1,0],
          [0,0,0,5.9]])
U=np.mat([[0,2.3,-4.4,0.5],
          [0,0,-3.5,2.8],
          [0,0,0,-1.3],
          [0,0,0,0]])
b=np.mat([[15.1],[1.8],[16.6],[36.9]])
def GS(a,b,L,D,U,n=10000,g=10e-7):
    if len(a) == len(b):  #若mx和mr长度相等则开始迭代 否则方程无解
        x = [] #迭代初值 初始化为单行全0矩阵
        c = []
        for i in range(len(b)):
            x.append([0])
            c.append([0])
    k=0
    ni=np.linalg.inv(D+L)
    print("中间运算所需逆阵:\n",ni)
    B=np.dot(ni,-U)
    f=np.dot(ni,b)
    print("中间运算矩阵B:\n",B)
    print("中间运算矩阵f:\n",f)
    while k<n:
        y=np.dot(B,x)+f
        c=np.mat(y-x)
        x=y
        zuida=c.max()
        print("c:\n",c)
        print("x:\n",x)
        print("误差c矩阵的最大值为：",c.max())
       
        if zuida<g:
            print("最终的x结果矩阵：\n",x)
            k=k+1
            print("循环次数：{cishu},运算完成".format(cishu=k))
            break
        else:
            k=k+1
            print("循环次数：{cishu}".format(cishu=k))
GS(a,b,L,D,U,n=10000,g=10e-7)