From 94c7e886bf5f84f73f5eb9fc0b863f8222db6138 Mon Sep 17 00:00:00 2001 From: lhye200 Date: Wed, 4 Jun 2025 21:39:19 +0800 Subject: [PATCH] 20250604 --- 227-8.py | 35 +++++++++++++++++++++ 228-12.py | 88 +++++++++++++++++++++++++++++++++++++++++++++++++++ 228-13.py | 56 +++++++++++++++++++++++++++++++++ 228-15.py | 60 +++++++++++++++++++++++++++++++++++ 228-17.py | 94 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 333 insertions(+) create mode 100644 227-8.py create mode 100644 228-12.py create mode 100644 228-13.py create mode 100644 228-15.py create mode 100644 228-17.py diff --git a/227-8.py b/227-8.py new file mode 100644 index 0000000..f72283b --- /dev/null +++ b/227-8.py @@ -0,0 +1,35 @@ + +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))] + + +if __name__ == "__main__": + A = [[1, 3], [-2, 4]] + 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=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)) + print("Norm of A with v=inf:", Norm(A, float("inf"))) + print("Norm of A with v=1:", Norm(A, 1)) + + \ No newline at end of file diff --git a/228-12.py b/228-12.py new file mode 100644 index 0000000..d09e170 --- /dev/null +++ b/228-12.py @@ -0,0 +1,88 @@ +def MullerSolve(fx,x0,x1,x2,err1,err2,N): + count = 0 + f0 = fx(x0) + f1 = fx(x1) + f2 = fx(x2) + q = (x2 - x1) / (x1 - x0) + p = 0 + a = 0 + b = 0 + c = 0 + while True: + p = (x2 - x0) / (x1 - x0) + + a = q**2 * f0 - q*p*f1 + q*f2 + b = q**2 *f0 - p**2 *f1 + (p + q)*f2 + c = p*f2 + h1 = 0 + + if b.real < 0: + h1 = -2 * c / (b - (b**2 - 4*a*c)**0.5) + else: + h1 = -2 * c / (b + (b**2 - 4*a*c)**0.5) + x3 = x2 + h1 * (x2 - x1) + f3 = fx(x3) + k = err1 + 1 + if abs(f3) < 1: + k = abs(x3 - x2) + else: + k = abs(x3 - x2) / abs(f3) + + if abs(f3) < err2 or k < err1: + return x3, 1 + count += 1 + if count > N: + return None, 0 + x0 = x1 + x1 = x2 + x2 = x3 + f0 = f1 + f1 = f2 + f2 = f3 + q = h1 + +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,0,1], + [2,2,1], + [-1,0,0] + ] + + lam = [] + + count = 0 + k = -100 + fx = lambda x: Det([[A[i][j] - x * (1 if i == j else 0) for j in range(len(A))] for i in range(len(A))]) + k = 10 + while len(lam) < len(A): + for i in range(-k,k): + re,sta = MullerSolve(fx, i, i + 1, i + 2, 1e-10, 1e-10, 100) + if sta == 1: + a = round(re.real, 9) + b = round(re.imag, 9) + re_t = complex(a, b) + if re_t not in lam: + if re_t.imag != 0: + lam.append(re_t) + lam.append(re_t.conjugate()) + else: + lam.append(a) + if len(lam) == len(A): + break + k *= 10 + p = abs(lam[0]) + for i in range(len(lam)): + print(f"λ{i+1} = {lam[i]}") + if abs(lam[i]) > p: + p = abs(lam[i]) + + print(f"谱半径 = {p:.3f}") \ No newline at end of file diff --git a/228-13.py b/228-13.py new file mode 100644 index 0000000..9e458b6 --- /dev/null +++ b/228-13.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 SOR(A,b,x,w,err,N): + count = 0 + n = len(A) + while True: + count += 1 + for i in range(n): + sum1 = sum(A[i][j] * x[j] for j in range(i)) + sum2 = sum(A[i][j] * x[j] for j in range(i, n)) + x[i] += w*(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__": + A = [ + [4,-1,0,-1,0,0], + [-1,4,-1,0,-1,0], + [0,-1,4,0,0,-1], + [-1,0,0,4,-1,0], + [0,-1,0,-1,4,-1], + [0,0,-1,0,-1,4] + ] + b = [2,3,2,2,1,2] + + x = [0,0, 0, 0, 0, 0] + err = 1e-5 + + w = 1 + x1,k,sta = SOR(A, b, x, w, err, 100) + print(f"w = {w}, 解为: {x1}, 迭代次数: {k}, 状态: {'收敛' if sta == 1 else '未收敛'}") + + w = 1.1 + x = [0, 0, 0, 0, 0, 0] + x2,k,sta = SOR(A, b, x, w, err, 100) + print(f"w = {w}, 解为: {x2}, 迭代次数: {k}, 状态: {'收敛' if sta == 1 else '未收敛'}") diff --git a/228-15.py b/228-15.py new file mode 100644 index 0000000..1693e0a --- /dev/null +++ b/228-15.py @@ -0,0 +1,60 @@ +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 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 + +def Inverse(A): + n = len(A) + # 计算代数余子式矩阵 + B = [[0 for i in range(n)] for j in range(n)] + for i in range(n): + for j in range(n): + minor = [row[:j] + row[j+1:] for row in (A[:i] + A[i+1:])] + B[j][i] = ((-1) ** (i + j)) * sum(minor[k][l] * (-1) ** (k + l) for k in range(n - 1) for l in range(n - 1)) + det = Det(A) + print(det) + if det == 0: + print("矩阵不可逆") + return None + A_inv = [[B[i][j] / det for j in range(n)] for i in range(n)] + return A_inv + +def Cond(A,v): + inv_A = Inverse(A) + print(inv_A,Norm(A, v), Norm(inv_A, v)) + return Norm(A, v) * Norm(inv_A, v) + + +if __name__ == "__main__": + A = [ + [1,2], + [1.001,2.001] + ] + + print(f"矩阵A的条件数为: {Cond(A, float('inf')):.5f}") + + A = [ + [1,2], + [3,4] + ] + print(f"矩阵A的条件数为: {Cond(A, float('inf')):.5f}") diff --git a/228-17.py b/228-17.py new file mode 100644 index 0000000..9399933 --- /dev/null +++ b/228-17.py @@ -0,0 +1,94 @@ +def SovleRowMain(A,b,round_num=15): + ks = 0.00000001 + n = len(A) + if len(A[0]) != n: + print("A要为方阵") + return None, None, None, None + if len(b) != n: + print("b与A的行数不匹配") + return None, None, None, None + p = list(range(n)) + for i in range(n): + row_max = abs(A[i][i]) + row_max_index = i + for j in range(i + 1, n): + if abs(A[j][i]) > row_max: + row_max = abs(A[j][i]) + row_max_index = j + A[i], A[row_max_index] = A[row_max_index], A[i] + b[i], b[row_max_index] = b[row_max_index], b[i] + p[i], p[row_max_index] = p[row_max_index], p[i] + + if abs(A[i][i]) < ks: + print("A矩阵奇异，无法进行高斯消元") + return None, None, None, None + for j in range(i + 1, n): + m = round(A[j][i] / A[i][i],round_num) + A[j][i] = m + for k in range(i + 1, n): + A[j][k] -= round(m * A[i][k],round_num) + b[j] -= round(m * b[i],round_num) + + if abs(A[n - 1][n - 1]) < ks: + print("A矩阵奇异，无法进行高斯消元") + return None, None, None, None + + # 回代求解 + b[n - 1] = round(b[n - 1]/A[n - 1][n - 1],round_num) + for i in range(n - 2, -1, -1): + for j in range(i + 1, n): + b[i] -= round(A[i][j] * b[j],round_num) + b[i] /= round(A[i][i]) + b = [round(b[i], round_num) for i in range(n)] + # 得到L,U和P矩阵 + L = [[0 for i in range(n)] for j in range(n)] + U = [[0 for i in range(n)] for j in range(n)] + P = [[0 for i in range(n)] for j in range(n)] + for i in range(n): + for j in range(n): + if i == j: + L[i][j] = 1 + U[i][j] = A[i][j] + elif i < j: + U[i][j] = A[i][j] + else: + L[i][j] = A[i][j] + + P[i][p[i]] = 1 + return P,L,U,b + +def IterativeMethod(A, b, err, N): + b_c = [b[i] for i in range(len(b))] + A_c = [[A[i][j] for j in range(len(A[0]))] for i in range(len(A))] + P,L,U,x0 = SovleRowMain(A_c, b_c,4) + print(L) + print(U) + print(f"初始解为: {x0}") + count = 0 + while count