合肥工业大学计算方法 Python 高分实验

import numpy as np def Gauss(A, b, x0, xstar, eps=0.0001, max_times=100): k = 0 while True: for i in range(0, len(A)): sum = 0.0 for j in range(0, i): sum = sum + A[i][j] * xstar[j] for j in range(i + 1, len(A)): sum = sum + A[i][j] * x0[j] xstar[i] = (b[i] - sum) / A[i][i] temp = np.fabs(xstar[0] - x0[0]) for j in range(1, len(A)): if np.fabs(xstar[j] - x0[j]) > temp: temp = np.fabs(xstar[j] - x0[j]) k = k + 1 if (temp < eps or k > max_times): print(k, end=" ") for i in range(len(xstar)): print("%.5f" % xstar[i], end=" ") print() print("误差限等于{0}时，高斯赛德尔迭代法需要迭代{1}次收敛".format(eps, k)) return None else: x0 = xstar.copy() print(k, end=" ") for i in range(len(xstar)): print("%.5f" % xstar[i], end=" ") print() def SOR(A, b, x0, xstar, w, eps=0.0001, max_times=100): # 超松弛法 k = 0 while True: for i in range(0, len(A)): sum = 0.0 for j in range(0, i): sum = sum + A[i][j] * xstar[j] for j in range(i + 1, len(A)): sum = sum + A[i][j] * x0[j] xstar[i] = (b[i] - sum) / A[i][i] xstar[i] = w * xstar[i] + (1 - w) * x0[i] temp = np.fabs(xstar[0] - x0[0]) for j in range(1, len(A)): if np.fabs(xstar[j] - x0[j]) > temp: temp = np.fabs(xstar[j] - x0[j]) k = k + 1 if (temp <= eps or k > max_times): print(k, end=" ") for i in range(len(xstar)): print("%.5f" % xstar[i], end=" ") print() print("误差限等于{0}时，超松弛法需要迭代{1}次收敛".format(eps, k)) return None else: x0 = xstar.copy() print(k, end=" ") for i in range(len(xstar)): print("%.5f" % xstar[i], end=" ") print() def Jacobi(A, b, x0, xstar, eps=0.0001, max_times=100): k = 0 while True: for i in range(0, len(A)): sum = 0.0 for j in range(0, i): sum = sum + A[i][j] * x0[j] for j in range(i + 1, len(A)): sum = sum + A[i][j] * x0[j] xstar[i] = (b[i] - sum) / A[i][i] temp = np.fabs(xstar[0] - x0[0]) for j in range(1, len(A)): if np.fabs(xstar[j] - x0[j]) > temp: temp = np.fabs(xstar[j] - x0[j]) k = k + 1 if (temp < eps or k > max_times): print(k, end=" ") for i in range(len(xstar)): print("%.5f" % xstar[i], end=" ") print() print("误差限等于{0}时，雅克比迭代法需要迭代{1}次收敛".format(eps, k)) return None else: x0 = xstar.copy() print(k, end=" ") for i in range(len(xstar)): print("%.5f" % xstar[i], end=" ") print() if __name__ == '__main__': # A = np.array([[10, -1, -2], [-1, 10, -2], [-1, -1, 5]]) # A= np.array([[3,1],[1,2]]) # A= np.array([[10,-1,2,0],[-1,11,-1,3],[2,-1,10,-1],[0,3,-1,8]]) # A= np.array([[3,1],[1,2]]) A = np.array([[4, 3, 0], [3, 4, -1], [0, -1, 4]]) # B = np.array([7.2, 8.3, 4.2]) # B =np.array([2,1]) # B =np.array([6,25,-11,15]) # B =np.array([-1,2]) B = np.array([16, 20, -12]) x0 = np.zeros([len(A), 1], dtype=float) x = np.zeros([len(A), 1], dtype=float) # x = np.array([0.0, 0, 0]) eps = 10 ** (-4) k = 100 # 最大迭代次数 print("k", end=" ") for i in range(len(A)): print("x{0}".format(i + 1), end=" ") print() Gauss(A, B, x0, x, eps, k) print("-------------------------") print("雅克比迭代") print("k", end=" ") for i in range(len(A)): print("x{0}".format(i + 1), end=" ") print() Jacobi(A, B, x0, x, eps, k) print("-------------------------") print("SOR") print("k", end=" ") for i in range(len(A)): print("x{0}".format(i + 1), end=" ") print() SOR(A, B, x0, x, 1.25,eps, k)