diff --git a/按方法整理/矩阵-模数-范数-点积.py b/按方法整理/矩阵-模数-范数-点积.py index a8f4fb6..c2a2567 100644 --- a/按方法整理/矩阵-模数-范数-点积.py +++ b/按方法整理/矩阵-模数-范数-点积.py @@ -25,7 +25,7 @@ if __name__ == "__main__": x = [[1], [-1]] # 注意储存形式,单独一列的相量,每个数字都要中括号 print("Norm of x with v=1:", Norm(x, 1)) - print("Norm of x with v=inf:", Norm(x, float("inf"))) + print("Norm of x with v=inf:", Norm(x, float("inf"))) # 1 1范数;2 2范数;float('inf') 无穷范数 print("Norm of x with v=2:", Norm(x, 2)) print("Dot product of A and x:", Dot(A, x)) print("Norm of Ax with v=2:", Norm(Dot(A, x), 2)) diff --git a/按方法整理/矩阵-雅可比得J.py b/按方法整理/矩阵-雅可比得J.py new file mode 100644 index 0000000..54af9a0 --- /dev/null +++ b/按方法整理/矩阵-雅可比得J.py @@ -0,0 +1,56 @@ +#模 范数 +def Norm(x,v): + if len(x[0]) == 1: + if v == 1: + return sum([abs(i[0]) for i in x]) + elif v == 2: + return (sum([i[0]**2 for i in x]))**0.5 + elif v == float("inf"): + return max([abs(i[0]) for i in x]) + else: + if v == 1: + return max([sum([abs(x[i][j]) for i in range(len(x))]) for j in range(len(x[0]))]) + elif v == float("inf"): + return max([sum([abs(i) for i in x[j]]) for j in range(len(x))]) + return None + +# 计算矩阵的点积 +def Dot(A,B): + if len(A[0]) != len(B): + return None + return [[sum([A[i][j] * B[j][k] for j in range(len(A[0]))]) for k in range(len(B[0]))] for i in range(len(A))] + +# 计算矩阵的行列式 +def Det(A): + if len(A) == 2: + return A[0][0] * A[1][1] - A[0][1] * A[1][0] + det = 0 + for c in range(len(A)): + sub_matrix = [row[:c] + row[c+1:] for row in A[1:]] + det += ((-1) ** c) * A[0][c] * Det(sub_matrix) + return det + + +if __name__ == "__main__": + ########################################################################## + # 把矩阵换成题干的矩阵 ######################### + A = [ + [1,2,-2], + [1,1,1], + [2,2,1] + ] + # 把范数的种类数换成题干的要求,inf是无穷范数 ######################### + LU = A.copy() + ID = [[0 for _ in range(len(A))] for __ in range(len(A))] + for i in range(len(A)): + for j in range(len(A[0])): + if i == j: + ID[i][j] = 1/A[i][j] + LU[i][j] = 0 + else: + LU[i][j] = -LU[i][j] + print(f"LU分解的LU矩阵为: {LU}") + print(f"LU分解的D矩阵的逆矩阵为: {ID}") + J = Dot(ID, LU) + print(f"LU分解的J矩阵为: {J}") + diff --git a/按方法整理/矩阵-雅可比迭代法.py b/按方法整理/矩阵-雅可比迭代法.py new file mode 100644 index 0000000..efc1b72 --- /dev/null +++ b/按方法整理/矩阵-雅可比迭代法.py @@ -0,0 +1,55 @@ +import math + +#模 范数 +def Norm(x,v): + if len(x[0]) == 1: + if v == 1: + return sum([abs(i[0]) for i in x]) + elif v == 2: + return (sum([i[0]**2 for i in x]))**0.5 + elif v == float("inf"): + return max([abs(i[0]) for i in x]) + else: + if v == 1: + return max([sum([abs(x[i][j]) for i in range(len(x))]) for j in range(len(x[0]))]) + elif v == float("inf"): + return max([sum([abs(i) for i in x[j]]) for j in range(len(x))]) + return None + +#Jacobi 雅可比迭代 +def Jacobi(A,b,x,err,N): + count = 0 + n = len(A) + while True: + count += 1 + x_tt = x.copy() # 保存上一次迭代的解 + for i in range(n): + sum1 = sum(A[i][j] * x_tt[j] for j in range(i)) + sum2 = sum(A[i][j] * x_tt[j] for j in range(i+1, n)) + x[i] = (b[i] - sum1 - sum2) / A[i][i] + r = [[b[i] - sum(A[i][j] * x[j] for j in range(len(A[0])))] for i in range(n)] + err_now = Norm(r, float("inf")) + x_t = [round(i,5) for i in x] + print(f"第{count}次迭代, 误差 = {err_now:.5}, x = {x_t}") + if err_now < err: + return x, count, 1 + + if count > N: + return None,count, 0 + + +if __name__ == "__main__": + ########################################################################## + #把矩阵改成题干的矩阵,b改成题干结果,err精度要求修改########################## + A = [ + [10,-1,2,0], + [-1,11,-1,3], + [2,-1,10,-1], + [0,3,-1,8] + ] + b = [6,25,-11,15] + + x = [0,0, 0, 0] # 初始解 + err = 1e-5 # 精度要求 + x1,k,sta = Jacobi(A, b, x, err, 100) + print(f"解为: {x1}, 迭代次数: {k}, 状态: {'收敛' if sta == 1 else '未收敛'}")