Merge branch 'main' of https://gitea.pi5.lhye.work/lhye200/CalWay_Python

120-3.py

@@ -4,6 +4,7 @@ x1, x2 = 3, 6
 def fx(x):
     return x/(4+x**2)
 
+# 复合Newton-Cotes公式 复合牛顿-科特斯公式
 # n等分参数，x1到x2的区间，type=1表示梯形法，type=2表示辛普森法
 def CompositeNewtonCotes(n, type):
     if type == 1:

120-4.py

@@ -2,6 +2,7 @@ def fx(x):
     return 1/x
 
 
+# 逐次分半梯形递推公式
 def SplitTrapezoidal(a,b,err):
     t1 = (b-a)*(fx(a)+fx(b))/2
     k = 1

120-5.py

@@ -5,6 +5,7 @@ def fx(x):
         x = 1e-10  # Avoid division by zero
     return math.sin(x)/x
 
+# 龙贝格方法 积分
 def Romberg(a, b, err):
     t00 = (b-a)*(fx(a)+fx(b))/2
     t01 = t10 = t11 = t20 = t21 = t30 = t31 = 0

159-1.py

@@ -1,8 +1,8 @@
 def fx(x):
     return x**4-3*x+1
 
-
-def erfen(x1,x2,err):
+# 二分法求解方程的根
+def SolveByDivTwo(x1,x2,err):
     while abs(x2-x1) >= err:
         x = (x1+x2)/2
         if fx(x) * fx(x1) < 0:
@@ -15,5 +15,6 @@ def erfen(x1,x2,err):
 if __name__ == "__main__":
     x1 = 0.3
     x2 = 0.4
-    root = erfen(x1, x2, 0.5e-5)
-    print(f"解为: {root}")
+    err = 0.5e-5
+    root = SolveByDivTwo(x1, x2, err)
+    print(f"Root: {root:.5f}")

159-2.py

@@ -4,14 +4,22 @@ import matplotlib.pyplot as plt
 def fx(x):
     return math.exp(x)-math.sin(x)
 
-def huatu(a, b, buchang):
-    x = [a + (b-a)*i*buchang for i in range(int(1/buchang+1))]
+# 绘制函数图像 并标记可能的零点
+def DrawGraph(a, b, stepper):
+    x = [a + (b-a)*i*stepper for i in range(int(1/stepper+1))]
     y = [fx(i) for i in x]
+    plt.xlim(a-abs(b-a)*0.1, b+abs(b-a)*0.1)
     plt.plot(x, y)
-    plt.axhline(0, color='black', lw=1)
+    plt.axhline(0, color='black', lw=0.5, ls='-')
+    plt.axvline(a, color='red', lw=0.5, ls='--')
+    plt.axvline(b, color='red', lw=0.5, ls='--')
+    plt.text(a, y[0], f'({a:.5f},{y[0]:.5f})', fontsize=8, ha='left')
+    plt.text(b, y[-1], f'({b:.5f},{y[-1]:.5f})', fontsize=8, ha='right')
     for i in range(len(x)-1):
         if y[i] * y[i+1] < 0:
-            print(f"({x[i]},{y[i]}),({x[i+1]},{y[i+1]})")
+            print(f"可能存在零点: ({x[i]:.5f},{y[i]:.5f})和({x[i+1]:.5f},{y[i+1]:.5f})之间")
+            plt.plot((x[i]+x[i+1])/2, (y[i]+y[i+1])/2, 'ro', markersize=3)
+    plt.title("Graph of f(x)")
     plt.xlabel("x")
     plt.ylabel("f(x)")
     plt.show()
@@ -21,4 +29,7 @@ def huatu(a, b, buchang):
 if __name__ == "__main__":
     a = -2*math.pi
     b = math.pi
-    huatu(a, b, 0.00001)
+    step = 0.00001
+    print(f"边界点: {a}, {b}")
+    print(f"步长: {step}")
+    x, y = DrawGraph(a, b, step)

159-3.py

@@ -1,40 +1,44 @@
-def f1(x):
+def fd1(x):
     return 1+1/x**2
 
-def f2(x):
-    return 1/((x-1)**0.5)
+def fd2(x):
+    return 1/(x-1)**0.5
 
-def diedai(x,fd,err):
+# 迭代法求解方程
+def Renew(x,fd,err):
     count=0
     i = 0
+    try:
         while True:
             xk = fd(x)
-        if complex(xk).imag != 0:
-            break
-        if (i) < abs(xk-x):
+            if abs(i) < (xk-x):
                 count += 1
             else:
                 count = 0
             if count > 10:
+                print("不收敛")
                 break
             if abs(xk-x) < err:
                 return xk
-        i = abs(xk-x)
+            i = xk-x
             x = xk
+    except Exception as e:
+        print(f"发生错误: {e}")
+        return None
     return None
 
 
 if __name__ == "__main__":
     x0 = 1.5
 
-    result = diedai(x0, f1, 1e-5)
+    result = Renew(x0, fd1, 1e-5)
     if result is not None:
-        print(f"(1)收敛 解为: {result:.5f}")
+        print(f"1式收敛 解为: {result:.5f}")
     else:
-        print("(1)不收敛")
+        print("1式不收敛")
 
-    result = diedai(x0, f2, 1e-5)
+    result = Renew(x0, fd2, 1e-5)
     if result is not None:
-        print(f"(2)收敛 解为: {result:.5f}")
+        print(f"2式收敛 解为: {result:.5f}")
     else:
-        print("(2)不收敛")
+        print("2式不收敛")

159-5.py

@@ -0,0 +1,61 @@
+import math
+
+
+def f1(x):
+    return x**2 + 10*math.cos(x)
+
+def df1(x):
+    return 2*x - 10*math.sin(x)
+
+
+def f2(x):
+    return 1 + math.atan(x) - x
+
+def df2(x):
+    return 1/(1+x**2) - 1
+
+# 牛顿方法求解方程
+def NewtonSolve(fx, dfx, x0, err, N0):
+    count = 0
+    x1 = x0 + 1 + err
+    while abs(x1 - x0) > err or abs(fx(x1)) > err: # 添加条件修正根误差太大的问题
+        if abs(dfx(x1)) < 1e-10:
+            return None, 0
+        x1 = x0 - fx(x0) / dfx(x0)
+        count += 1
+        if count > N0:
+            return None, -1
+        x0 = x1
+    return x1, 1
+
+# 查找根区间
+def FindRootZone(fx,start,stop,step):
+    x = start
+    while x < stop:
+        if fx(x) * fx(x+step) < 0:
+            return x
+        x += step
+    return None
+
+if __name__ == "__main__":
+    err = 1e-5
+    N0 = 100
+
+    x0 = 1.6
+    result,status = NewtonSolve(f1, df1, x0, err,N0)
+    if status == 1:
+        print(f"f1收敛 解为: {result}")
+    elif status == -1:
+        print("f1不收敛")
+    else:
+        print("f1导数为0，无法收敛")
+
+
+    x0 = FindRootZone(f2, -5, 5, 0.01)
+    result, status = NewtonSolve(f2, df2, x0, err,N0)
+    if status == 1:
+        print(f"f2收敛 解为: {result}")
+    elif status == -1:
+        print("f2不收敛")
+    else:
+        print("f2导数为0，无法收敛")

159-6.py

@@ -0,0 +1,39 @@
+import math
+
+
+def f(x):
+    return x**2 - 30
+
+def df(x):
+    return 2*x
+
+# 牛顿方法求解方程
+def NewtonSolve(fx, dfx, x0, err, N0):
+    count = 0
+    x1 = x0 + 1 + err
+    while abs(x1 - x0) > err or abs(fx(x1)) > err: # 添加条件修正根误差太大的问题
+        if abs(dfx(x1)) < 1e-10:
+            return None, 0
+        x1 = x0 - fx(x0) / dfx(x0)
+        count += 1
+        if count > N0:
+            return None, -1
+        x0 = x1
+    return x1, 1
+
+# 查找根区间
+def FindRootZone(fx,start,stop,step):
+    x = start
+    while x < stop:
+        if fx(x) * fx(x+step) < 0:
+            return x
+        x += step
+    return None
+
+if __name__ == "__main__":
+    err = 1e-4
+    N0 = 100
+
+    x0 = FindRootZone(f, 5, 6, 0.01)
+    result,status = NewtonSolve(f, df, x0, err,N0)
+    print(f"sqrt(30) = {result:.3f}")

159-7.py

@@ -0,0 +1,37 @@
+# 查找根区间
+def FindRootZone(fx,start,stop,step):
+    x = start
+    while x < stop:
+        if fx(x) * fx(x+step) < 0:
+            return x
+        x += step
+    return None
+
+# 求n次方根迭代过程如下
+def GetNthRoot(a,n):
+    if a < 0 and n % 2 == 0:
+        print("Cannot compute even root of negative number")
+        return None
+    fx = lambda x: x**n - a
+    dfx = lambda x: n * x**(n-1)
+    err = 1e-10
+    N0 = 100
+    x0 = 0
+    if a > 0:
+        x0 = FindRootZone(fx, 0, a, 0.01)
+    else:
+        x0 = FindRootZone(fx, a, 0, 0.01)
+
+    count = 0
+    x1 = x0 + 1 + err
+    while abs(x1 - x0) > err or abs(fx(x1)) > err: # 添加条件修正根误差太大的问题
+        x1 = x0 - fx(x0) / dfx(x0)
+        count += 1
+        if count > N0:
+            return None
+        x0 = x1
+    return x1
+
+if __name__ == "__main__":
+    re = GetNthRoot(30, 5)
+    print(re)

160-11.py

@@ -0,0 +1,35 @@
+#正割法计算方程
+def SecantSolve(fx, x0, x1, err=1e-10, N0=100):
+    count = 0
+    print(f"k={count}: x{count}={x0}, f(x{count})={fx(x0)}")
+    count += 1
+    print(f"k={count}: x{count}={x1}, f(x{count})={fx(x1)}")
+    count += 1
+    while abs(x1 - x0) > err or abs(fx(x1)) > err:
+        if fx(x1) == fx(x0):
+            return None,0
+        x2 = x1 - fx(x1) * (x1 - x0) / (fx(x1) - fx(x0))
+        print(f"k={count}: x{count}={x2}, f(x{count})={fx(x2)}")
+        count += 1
+        if count > N0:
+            return None,-1
+        x0 = x1
+        x1 = x2
+    return x2,1
+
+
+if __name__ == "__main__":
+    err = 1e-5
+    N0 = 100
+
+    x0 = 0.3
+    x1 = 0.4
+    fx = lambda x: x**4 - 3*x + 1
+
+    result, status = SecantSolve(fx, x0, x1, err, N0)
+    if status == 1:
+        print(f"fx收敛 解为: {result}")
+    elif status == -1:
+        print("fx不收敛")
+    else:
+        print("分母为0，无法收敛")

160-12.py

@@ -0,0 +1,63 @@
+#抛物线法解方程（主）
+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 < 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)
+
+        print(f"k={count}: x{count}={x0:.5}, f(x{count})={f0:.5}; x{count+1}={x1:.5}, f(x{count+1})={f1:.5}; x{count+2}={x2:.5}, f(x{count+2})={f2:.5}; x{count+3}={x3:.5}, f(x{count+3})={f3:.5}")
+        print(f"p={p:.5}, q={q:.5}, a={a:.5}, b={b:.5}, c={c:.5}, h={h1:.5}")
+
+        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
+
+
+if __name__ == "__main__":
+    err1 = 1e-5
+    err2 = 1e-5
+    N = 100
+
+    x0 = 0.3
+    x1 = 0.5
+    x2 = 0.4
+    fx = lambda x: 8*x**4 - 8*x**2 + 1
+
+    result, status = MullerSolve(fx, x0, x1, x2, err1, err2, N)
+    if status == 1:
+        print(f"fx收敛 解为: {result}")
+    else:
+        print("fx不收敛")

227-3.py

@@ -1,10 +1,3 @@
-A = [
-    [0, 3, 4],
-    [1, -1, 1],
-    [2, 1, 2]
-]
-
-b = [1, 2, 3]
 
 # 列主元高斯消元法
 def SovleRowMain(A,b):
@@ -61,7 +54,7 @@ def SovleRowMain(A,b):
 
         P[i][p[i]] = 1
     return P,L,U,b
-
+#打印矩阵
 def prettyPrintMatrix(matrix):
     for row in matrix:
         print(row)
@@ -70,6 +63,15 @@ def prettyPrintMatrix(matrix):
 if __name__ == "__main__":
     # A = np.array(A, dtype=float)
     # b = np.array(b, dtype=float)
+      
+    A = [
+        [0, 3, 4],
+        [1, -1, 1],
+        [2, 1, 2]
+    ]
+
+    b = [1, 2, 3]
+
     P,L,U,x = SovleRowMain(A, b)
     print("P:")
     prettyPrintMatrix(P)
@@ -80,3 +82,4 @@ if __name__ == "__main__":
     print("x:")
     print(x)
         
+ 

227-4.py

@@ -1,8 +1,3 @@
-# 储存下三角矩阵 a11, a21, a22, a31, a32, a33 ...
-A = [4,2,2,-2,-3,14]
-
-
-b = [10,5,4]
 
 #列 行
 def getIndexFromDownMatrix(col, row):
@@ -41,5 +36,11 @@ def SqrtSolve(A,b):
 
 
 if __name__ == "__main__":
+    # 储存下三角矩阵 a11, a21, a22, a31, a32, a33 ...
+    A = [4,2,2,-2,-3,14]
+
+
+    b = [10,5,4]
+
     print("x:")
     print(SqrtSolve(A,b))

227-7.py

@@ -1,14 +1,3 @@
-# 储存追赶法A矩阵
-A = [
-    [0,4,-1],
-    [-1,4,-1],
-    [-1,4,-1],
-    [-1,4,-1],
-    [-1,4,0]
-]
-
-
-b = [100,200,200,200,100]
 
 # 追赶法求解
 def ZGsolve(A,b):
@@ -33,5 +22,16 @@ def ZGsolve(A,b):
 
 
 if __name__ == "__main__":
+    # 储存追赶法A矩阵
+    A = [
+        [0,4,-1],
+        [-1,4,-1],
+        [-1,4,-1],
+        [-1,4,-1],
+        [-1,4,0]
+    ]
+
+
+    b = [100,200,200,200,100]
     print("x:")
     print(ZGsolve(A,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))
+    
+    

228-12.py

@@ -0,0 +1,89 @@
+#抛物线法解方程
+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}")

228-13.py

@@ -0,0 +1,57 @@
+#模 范数
+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
+#SOR方法 逐次超松弛迭代
+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 '未收敛'}")

228-15.py

@@ -0,0 +1,61 @@
+#模 范数
+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}")

228-17.py

@@ -0,0 +1,95 @@
+#列主元高斯消元法
+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<N:
+        r1 = [b[i] - sum([A[i][j] * x0[j] for j in range(len(A[0]))]) for i in range(len(A))]
+        A_c = [[A[i][j] for j in range(len(A[0]))] for i in range(len(A))]
+        d1 = SovleRowMain(A_c, r1,4)[3]
+        x0 = [x0[i] + d1[i] for i in range(len(x0))]
+        print(f"第{count+1}次迭代, x{count+2} = {x0}, r{count+1} = {r1}, d{count+1} = {d1}")
+        err_now = max(abs(r1[i]) for i in range(len(r1)))
+        count += 1
+        if err_now < err:
+            break
+    return x0,count
+
+
+if __name__ == "__main__":
+    A = [
+        [51,82],
+        [151/3,81]
+    ]
+    b = [235,232]
+
+    err = 1e-4
+    N = 1000
+
+    x = IterativeMethod(A, b, err, N)[0]
+    print(f"解为: {x}")
+
+

279-3.py

@@ -0,0 +1,26 @@
+#经典R-K，龙格-库塔法
+def ClassicRK(x0,y0,h,xk,fxy):
+    k1=k2=k3=k4=0
+    result = [(x0,y0)]
+    while x0<=xk:
+        k1 = fxy(x0,y0)
+        k2 = fxy(x0+h/2,y0+h*k1/2)
+        k3 = fxy(x0+h/2,y0+h*k2/2)
+        k4 = fxy(x0+h,y0+h*k3)
+        y0 += h*(k1+4*k2+k3)/6
+        x0 += h
+        result.append((x0,y0))
+    return result
+
+
+if __name__=="__main__":
+    x0 = 0
+    y0 = -1
+    fxy = lambda x,y: x + y
+    real_fx = lambda x: -x-1
+    h = 0.1
+    xk = 2
+    result = ClassicRK(x0, y0, h, xk, fxy)
+    print("x\ty\t\treal_y")
+    for x, y in result:
+        print(f"{x:.2f}\t{y}\t{real_fx(x)}")

68-1.py

@@ -1,3 +1,4 @@
+
 import math
 
 list_x = [10,11,12,13]

69-5.py

@@ -1,6 +1,7 @@
 x = [0,1,2,3]
 y = [0,0,0,0]
 
+# 追赶法
 def ZGsolve(A,b):
     n = len(b)
     beta = [0]*n
@@ -21,19 +22,21 @@ def ZGsolve(A,b):
 
     return b
 
-
+# 获取相邻点的差分
 def GetDList(list_r):
     result = []
     for i in range(1,len(list_r)):
         result.append(list_r[i] - list_r[i-1])
     return result
 
+# 获取相邻点的差商
 def GetDQList(list_x, list_y):
     result = []
     for i in range(1,len(list_y)):
         result.append((list_y[i] - list_y[i-1]) / (list_x[i] - list_x[i-1]))
     return result
 
+# 三次样条插值
 def CubicSplineInterpolation(list_x,list_y,boundary_type,a1,a2):
     list_h = GetDList(list_x)
     list_dqxy = GetDQList(list_x, list_y)
@@ -68,6 +71,7 @@ def CubicSplineInterpolation(list_x,list_y,boundary_type,a1,a2):
 
     return M,list_h
 
+# 打印矩阵
 def PrintResult(M,list_h,list_x,list_y):
     for i in range(len(list_h)):
         k1 = (M[i+1]-M[i])/6/list_h[i]

89-1.py

@@ -43,7 +43,7 @@ def SovleRowMain(A,b):
     return b
 
 
-
+# 最小二乘法拟合
 def LeastSquares(list_x,list_y,n):
     m = len(list_x)
     x_n = []
@@ -63,6 +63,7 @@ def LeastSquares(list_x,list_y,n):
         A.append(tmp)
     return SovleRowMain(A,b)
 
+# 计算多项式在给定x值上的值
 def CalculateY(list_x, coeff):
     re = []
     for i in range(len(list_x)):
@@ -71,6 +72,7 @@ def CalculateY(list_x, coeff):
             re[i] += coeff[j]*list_x[i]**j
     return re
 
+# 计算均方根误差
 def MeanSquareErr(list_y,list_y_approx):
     m = len(list_y)
     err = 0
@@ -78,6 +80,7 @@ def MeanSquareErr(list_y,list_y_approx):
         err += (list_y[i] - list_y_approx[i])**2
     return err**0.5
 
+# 计算最大误差
 def MaxErr(list_y,list_y_approx):
     m = len(list_y)
     err = 0
@@ -86,6 +89,7 @@ def MaxErr(list_y,list_y_approx):
             err = abs(list_y[i] - list_y_approx[i])
     return err
 
+# 打印拟合方程
 def PrintEquation(coeff):
     n = len(coeff)
     str_ = str(coeff[0]) + "+"

89-2.py

@@ -45,7 +45,7 @@ def SovleRowMain(A,b):
     return b
 
 
-
+# 最小二乘法拟合
 def LeastSquares(list_x,list_y,n):
     m = len(list_x)
     x_n = []

89-3.py

@@ -43,7 +43,7 @@ def SovleRowMain(A,b):
     return b
 
 
-
+# 最小二乘法拟合
 def LeastSquares(list_x,list_y,n):
     m = len(list_x)
     x_n = []