ee

按方法整理/矩阵-模数-范数-点积.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))

按方法整理/矩阵-雅可比得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}")
+    

按方法整理/矩阵-雅可比迭代法.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 '未收敛'}")