解线性方程组的迭代法：Jacobi迭代法与Gauss-Seide

Jacobi迭代法：

import numpy as np

# Jacobi 迭代法计算线性方程

# Ax = b, x_0为初始值， epsilon为误差限， n为最大迭代次数

A = np.array([[1,2,-2],

            [1,1,1],

            [2,2,1]])

b = np.array([1,2,3])

x_0 = np.array([0,0,0])

epsilon = 0.01

n = 100

def jacobi(A, b, x_0, n, epsilon):

    s = len(A)

    D = np.diag(np.diag(A))

    inv = np.linalg.inv(D)

    M = np.eye(s)-np.dot(inv,A)

    count = 0

    # rho 为谱半径

    rho = np.max(np.linalg.eig(M)[0])

    if rho >= 1:

        print("Jacobi迭代法发散")

    else:

        x = x_0

        err = 0x3f3f3f

        while 1:

            if count >= n:

                break

            if err < epsilon:

                break

            x = np.dot(M, x_0)+np.dot(inv, b)

            print("第",count,"次迭代结果：",x[0],',',x[1],',',x[2],'\n')

            count = count + 1

            err = abs(max(x-x_0))

            x_0 = ximport numpy as np

# Gauss-Seidel 迭代法计算线性方程

# Ax = b, x_0为初始值， epsilon为误差限， n为最大迭代次数

A = np.array([[1,2,-2],

            [1,1,1],

            [2,2,1]])

b = np.array([1,2,3])

x_0 = np.array([0,0,0])

epsilon = 0.01

n = 100

def GS(A, b, x_0, n, epsilon):

    D = np.diag(np.diag(A))

    L = -1*np.tril(A-D)

    U = -1*np.triu(A-D)

    inv = np.linalg.inv(D-L)

    M = np.dot(inv,U)

    count = 0

    # rho 为谱半径

    rho = abs(np.max(np.linalg.eig(M)[0]))

    if rho >= 1:

        print("Gauss-Seidel迭代法发散")

    else:

        x = x_0

        err = 0x3f3f3f

        while 1:

            if count >= n:

                break

            if err < epsilon:

                break

            x = np.dot(M, x_0)+np.dot(inv, b)

            print("第",count,"次迭代结果：",x[0],',',x[1],',',x[2],'\n')

            count = count + 1

            err = abs(max(x-x_0))

            x_0 = x

    print("谱半径:",rho)

GS(A, b.T, x_0, n, epsilon)

        print("谱半径:",rho)

jacobi(A, b.T, x_0, n, epsilon)

Gauss-Seidel迭代法：