「深度学习福利」大神带你进阶工程师，立即查看>>>
定义矩阵乘法 def mult(h, x): result = [] for x in h: summ = 0 for index, y in enumerate(x): summ += y * x[index] result.append(summ) return result
创建希尔伯特矩阵 def create_hobert(n): h=[] for x in range(1, n + 1): row = [] for y in range(1, n + 1): row.append(1 / (x + y - 1)) h.append(row) return h
雅克比迭代法 def jacobi(A,B,sigma,N): n = len(A) x0 = [] x = [] for i in range(0,n): x0.append(0) x.append(0) for k in range(1,N+1): R = 0 for i in range(0,n): sum_ax = 0 for j in range(0,n): sum_ax = sum_ax + A[i][j] * x0[j] x[i] = x0[i] + ((B[i] - sum_ax)/A[i][i]) if abs(x[i] - x0[i]) > R: R = abs(x[i] - x0[i]) if R <= sigma: print("精确度等于",sigma,"时，jacobi迭代法需要迭代",k,"次收敛") return (x,k) for i in range(0,n): x0[i] = x[i] return (x,k)
Hauss-Seidel迭代法 def gauss_seidel(A,B,sigma,N): n = len(A) x0 = [] x = [] for i in range(0,n): x0.append(0) x.append(0) for k in range(1,N+1): R = 0 for i in range(0,n): sum_ax = 0 for j in range(0,n): if j >= i: sum_ax = sum_ax + A[i][j] * x0[j] else: sum_ax = sum_ax + A[i][j] * x[j] x[i] = x0[i] + ((B[i] - sum_ax)/A[i][i]) if abs(x[i] - x0[i]) > R: R = abs(x[i] - x0[i]) if R <= sigma: print("精确度等于",sigma,"时，gauss_seidel迭代法需要迭代",k,"次收敛") return (x,k) for i in range(0,n): x0[i] = x[i] return (x,k)
SOR迭代法 def sor(A,B,sigma,N,omega): n = len(A) x0 = [] x = [] for i in range(0,n): x0.append(0) x.append(0) for k in range(1,N+1): R = 0 for i in range(0,n): sum_ax = 0 for j in range(0,n): if j >= i: sum_ax = sum_ax + A[i][j] * x0[j] else: sum_ax = sum_ax + A[i][j] * x[j] x[i] = x0[i] + omega * ((B[i] - sum_ax)/A[i][i]) if abs(x[i] - x0[i]) > R: R = abs(x[i] - x0[i]) if R <= sigma: print("精确度等于",sigma,"时，松弛因子为",omega,"时,sor迭代法需要迭代",k,"次收敛") print(x) return (x,k) for i in range(0,n): x0[i] = x[i] return (x,k)
测试希尔伯特矩阵H if __name__ == "__main__": n = 10 H = create_hobert(n) x = [1 for x in range(n)] b = mult(H,x) print('雅克比迭代法:',jacobi(H, b, 0.1, 200)) print('SOR迭代法:',sor(H, b, 0.001, 20000, 1.5)) # print('Guass-Seidel迭代法:',gauss_seidel(H,b,0.00001,20))