diff --git a/numerical_analysis/1/README.md b/numerical_analysis/1/README.md index d493de7..90ff896 100644 --- a/numerical_analysis/1/README.md +++ b/numerical_analysis/1/README.md @@ -56,5 +56,8 @@ end 作业提交有效时间是今天到10月23日(两周后)之前的任意时间。提交作业请将代码和报告打包,以“课后作业1-名字-学号”命名提交。 +![image-20210116232529898](https://public.veypi.com/img/screenshot/20210116232529.png) +![image-20210116235118221](https://public.veypi.com/img/screenshot/20210116235118.png) +![image-20210116235324496](https://public.veypi.com/img/screenshot/20210116235324.png) \ No newline at end of file diff --git a/numerical_analysis/1/main.py b/numerical_analysis/1/main.py index ac2ddfe..b9b9915 100644 --- a/numerical_analysis/1/main.py +++ b/numerical_analysis/1/main.py @@ -28,7 +28,7 @@ def LU_decomposition(A): return L, U -# 生成范围在 [-100, 100]之前的方阵A 和 b +# 生成随机矩阵 A, b, 范围[-10, 10] def randomAb(m): A = np.random.random([m, m]) * 20 - 10 dia = np.random.random(m) * 10 @@ -37,6 +37,29 @@ def randomAb(m): return A, np.random.randint(0, 10, [m, 1]) +# 生成稀疏矩阵A, b +def sparseMatrix(m): + A, b = randomAb(m) + for i in range(m): + for j in range(m): + if i != j: + a = A[i][j] + if abs(a) < 9: + A[i][j] = 0 + elif a > 0: + A[i][j] = 10 * (a - 9) + else: + A[i][j] = 10 * (9 + a) + + return A, b + + +# 生成病态矩阵A, b +def illMatrix(m): + A, b = randomAb(m) + return A, b + + class Question1: """ 求解 Ax=b @@ -104,60 +127,58 @@ class Question1: """ return np.linalg.solve(A, b) - def RMSE(self, solver, n=8): + def RMSE(self, solver, n=8, randFunc=randomAb): ## 计算方差 s = time.time() - A, b = randomAb(n) + A, b = randFunc(n) X = (A.dot(solver(A, b)) - b) ** 2 for i in range(1000): - A, b = randomAb(n) + A, b = randFunc(n) X = X + (A.dot(solver(A, b)) - b) ** 2 return np.max(X), time.time() - s def show(self): - N = [2 ** i for i in range(12)] - Y = [[0 for i in range(12)] for _ in range(4)] - Z = [[0 for i in range(12)] for _ in range(4)] + n = 12 + N = [2 ** i for i in range(n)] + Y = [[0 for i in range(n)] for _ in range(4)] + Z = [[0 for i in range(n)] for _ in range(4)] - plt.subplot(1, 2, 1) - for i in range(len(N)): - print("size: %s" % N[i]) - print("LU") - Y[0][i], Z[0][i] = self.RMSE(self.solver1, N[i]) - print("jacobi") - Y[1][i], Z[1][i] = self.RMSE(self.solver3, N[i]) - print("inverse") - Y[2][i], Z[2][i] = self.RMSE(self.solver3, N[i]) - print("default\n") - Y[3][i], Z[3][i] = self.RMSE(self.solver4, N[i]) + randomFuns = [randomAb, sparseMatrix, illMatrix] + for r in range(3): + plt.subplot(3, 2, 2*r + 1) + for i in range(n): + print("size: %s" % N[i]) + print("LU") + Y[0][i], Z[0][i] = self.RMSE(self.solver1, N[i]) + print("jacobi") + Y[1][i], Z[1][i] = self.RMSE(self.solver3, N[i]) + print("inverse") + Y[2][i], Z[2][i] = self.RMSE(self.solver3, N[i]) + print("default\n") + Y[3][i], Z[3][i] = self.RMSE(self.solver4, N[i]) - # N = range(12) - plt.plot(N, Y[0], label="LU") - plt.plot(N, Y[1], label="Jacobi") - plt.plot(N, Y[2], label="inverse") - plt.plot(N, Y[3], label="default solver") - # plt.xticks([0, 10, 100, 1000], [0, 10, 100, 1000]) - plt.title('Accuracy') - plt.yscale('symlog') - plt.xscale('symlog') - plt.legend(loc='lower right') - plt.subplot(1, 2, 2) - plt.plot(N, Z[0], label="LU") - plt.plot(N, Z[1], label="Jacobi") - plt.plot(N, Z[2], label="inverse") - plt.plot(N, Z[3], label="default solver") - plt.xscale('symlog') - plt.yscale('symlog') - plt.title('time cost') - plt.legend(loc='lower right') + # N = range(12) + plt.plot(N, Y[0], label="LU") + plt.plot(N, Y[1], label="Jacobi") + plt.plot(N, Y[2], label="inverse") + plt.plot(N, Y[3], label="default solver") + # plt.xticks([0, 10, 100, 1000], [0, 10, 100, 1000]) + plt.title('Accuracy') + plt.yscale('symlog') + plt.xscale('symlog') + # plt.legend(loc='lower right') + plt.subplot(3, 2, 2 * r + 2) + plt.plot(N, Z[0], label="LU") + plt.plot(N, Z[1], label="Jacobi") + plt.plot(N, Z[2], label="inverse") + plt.plot(N, Z[3], label="default solver") + plt.xscale('symlog') + plt.yscale('symlog') + plt.title('time cost') + # plt.legend(loc='lower right') plt.show() if __name__ == "__main__": q = Question1() - A, b = randomAb(3) - print(A) - print(np.linalg.det(A)) - # print(q.solver2(A, b)) - # print(q.solver4(A, b)) q.show()