Jacobi迭代法：
import numpy as np
# Jacobi 迭代法計算線性方程
# Ax = b, x_0為初始值， epsilon為誤差限， n為最大迭代次數(shù)
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迭代法發(fā)散")
? ? 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為最大迭代次數(shù)
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迭代法發(fā)散")
? ? 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迭代法：
