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迭代法：
网友评论