df

120-3.py

@@ -1,10 +1,11 @@
+#把范围与原函数换成题干的############
 x1, x2 = 3, 6
 
-
 def fx(x):
     return x/(4+x**2)
 
-# n等分参数，x1到x2的区间，type=1表示梯形法，type=2表示辛普森法
+# 复合Newton-Cotes公式 复合牛顿-柯特斯公式
+# n等分参数，x1到x2的区间，type=1表示梯形法，type=2表示辛普生法
 def CompositeNewtonCotes(n, type):
     if type == 1:
         h = (x2 - x1) / n
@@ -24,8 +25,8 @@ def CompositeNewtonCotes(n, type):
 
 if __name__ == "__main__":
     # 复合梯形公式，点数为n+1
-    print("复合梯形公式\n", CompositeNewtonCotes(8, 1))
+    print("复合梯形公式\n", CompositeNewtonCotes(8, 1)) #8等分，1代表是梯形公式####################
-    # 复合辛普森公式，点数为2n+1
+    # 复合辛普生公式，点数为2n+1
-    print("复合辛普森公式\n", CompositeNewtonCotes(4, 2))
+    print("复合辛普生公式\n", CompositeNewtonCotes(4, 2)) #4等分，2代表是辛普生公式###############
 
 

120-4.py

@@ -1,22 +1,28 @@
+#把函数改成题干的####################
 def fx(x):
     return 1/x
 
 
+# 逐次分半梯形递推公式
 def SplitTrapezoidal(a,b,err):
+    count = 1
     t1 = (b-a)*(fx(a)+fx(b))/2
+    print(f"t{count}={t1}")
     k = 1
     while True:
         tmp = 0
         for i in range(1, 2**(k-1)+1):
             tmp += fx(a + (b-a)*(2*i-1)/(2**k))
         t2 = t1/2+(b-a)*tmp/(2**k)
+        count *= 2
+        print(f"t{count}={t2}")
         if abs(t2-t1) < err:
             break
         t1 = t2
         k += 1
     return t2,k
 
-
+#把范围ab换成题干的范围，err换成题干的精度要求############
 if __name__ == "__main__":
     a = 1
     b = 3

120-5.py

@@ -2,11 +2,16 @@ import math
 
 def fx(x):
     if x == 0:
-        x = 1e-10  # Avoid division by zero
+        x = 1e-10  # Avoid division by zero  #如果x能为0，注释掉这行##############
-    return math.sin(x)/x
+        pass
+    return math.sin(x)/x   #把函数改成题干的形式###################
 
+# 龙贝格方法 积分
 def Romberg(a, b, err):
+    table = [[],[],[],[]]
     t00 = (b-a)*(fx(a)+fx(b))/2
+    table[0].append(t00)
+    print(f"t0(0)={t00}")
     t01 = t10 = t11 = t20 = t21 = t30 = t31 = 0
     k = 1
     while True:
@@ -14,34 +19,49 @@ def Romberg(a, b, err):
         for i in range(1, 2**(k-1)+1):
             tmp += fx(a + (b-a)*(2*i-1)/(2**k))
         t01 = t00/2+(b-a)*tmp/(2**k)
+        print(f"t0({k})={t01}")
+        table[0].append(t01)
         
         if k>1:
             t11 = (4*t01-t00)/3
-            
+            print(f"t1({k-1})={t11}")
+            table[1].append(t11)
             if k>2:
                 t21 = (16*t11-t10)/15
-                
+                print(f"t2({k-2})={t21}")
+                table[2].append(t21)
                 if k>3:
                     t31 = (64*t21-t20)/63
+                    print(f"t3({k-3})={t31}")
+                    table[3].append(t31)
                     if abs(t31-t30) < err:
                         break
                     t30 = t31
                 else:
                     t30 = (64*t21-t20)/63
+                    print(f"t3(0)={t30}")
+                    table[3].append(t30)
                 t20 = t21
             else:
                 t20 = (16*t11-t10)/15
+                print(f"t2(0)={t20}")
+                table[2].append(t20)
             t10 = t11
         else:
             t10 = (4*t01-t00)/3
+            print(f"t1(0)={t10}")
+            table[1].append(t10)
         t00 = t01
         k += 1
+    print("Romberg table:")
+    for i in range(len(table)):
+        print(f"t{i}: {table[i]}")
     return t31, k
 
-
+# 把范围ab与精度换成题干的############
 if __name__ == "__main__":
     a = 0
     b = 1
-    err = 0.5e-6
+    err = 0.5e-6    #把精度要求改成题干的##################
     result, k = Romberg(a, b, err)
     print(f"Result: {result}, k={k}")

159-1.py

@@ -1,20 +1,26 @@
+#原函数换成题干的*********************
 def fx(x):
     return x**4-3*x+1
 
-
+# 二分法求解方程的根
 def SolveByDivTwo(x1,x2,err):
+    count = 1
+    
     while abs(x2-x1) >= err:
         x = (x1+x2)/2
+        print(f"k={count},a{count}={x1},b{count}={x2},x{count}={x},fx(x)={fx(x)}")
         if fx(x) * fx(x1) < 0:
             x2 = x
         else:
             x1 = x
+        count += 1
+
     return (x1+x2)/2
 
-
+#范围和精度要求换成题干的############
 if __name__ == "__main__":
     x1 = 0.3
     x2 = 0.4
-    err = 0.5e-5
+    err = 0.5e-5     #精度要求##################
     root = SolveByDivTwo(x1, x2, err)
     print(f"Root: {root:.5f}")  

159-2.py

@@ -1,9 +1,11 @@
 import math
 import matplotlib.pyplot as plt
 
+#把原函数换成题干的形式########################
 def fx(x):
     return math.exp(x)-math.sin(x)
 
+# 绘制函数图像 并标记可能的零点
 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]
@@ -24,7 +26,7 @@ def DrawGraph(a, b, stepper):
     plt.show()
     return x, y
 
-
+#把范围改成题干的形式#######################
 if __name__ == "__main__":
     a = -2*math.pi
     b = math.pi

159-3.py

@@ -4,6 +4,7 @@ def fd1(x):
 def fd2(x):
     return 1/(x-1)**0.5
 
+# 迭代法求解方程
 def Renew(x,fd,err):
     count=0
     i = 0

159-5.py

@@ -1,6 +1,6 @@
 import math
 
-
+#原函数和导数改成题干的形式#####################
 def f1(x):
     return x**2 + 10*math.cos(x)
 
@@ -14,20 +14,23 @@ def f2(x):
 def df2(x):
     return 1/(1+x**2) - 1
 
-
+# 牛顿方法求解方程
 def NewtonSolve(fx, dfx, x0, err, N0):
     count = 0
+    print(f"k={count}, x0={x0}")
     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
+        print(f"k={count}, x{count}={x1},x1-x0={abs(x1-x0)}")
         if count > N0:
             return None, -1
         x0 = x1
     return x1, 1
 
+# 查找根区间
 def FindRootZone(fx,start,stop,step):
     x = start
     while x < stop:
@@ -40,6 +43,7 @@ if __name__ == "__main__":
     err = 1e-5
     N0 = 100
 
+#把初始值换成题干形式###############
     x0 = 1.6
     result,status = NewtonSolve(f1, df1, x0, err,N0)
     if status == 1:
@@ -50,7 +54,7 @@ if __name__ == "__main__":
         print("f1导数为0，无法收敛")
 
 
-    x0 = FindRootZone(f2, -5, 5, 0.01)
+    x0 = FindRootZone(f2, -5, 5, 0.01) # 查找f2的根区间与步长
     result, status = NewtonSolve(f2, df2, x0, err,N0)
     if status == 1:
         print(f"f2收敛 解为: {result}")

159-6.py

@@ -1,26 +1,29 @@
 import math
 
-
+#把函数改成题干的形式#####################
 def f(x):
     return x**2 - 30
 
 def df(x):
     return 2*x
 
-
+# 牛顿方法求解方程
 def NewtonSolve(fx, dfx, x0, err, N0):
     count = 0
+    print(f"k={count}, x0={x0}")
     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
+        print(f"k={count}, x{count}={x1},x1-x0={abs(x1-x0)}")
         if count > N0:
             return None, -1
         x0 = x1
     return x1, 1
 
+# 查找根区间
 def FindRootZone(fx,start,stop,step):
     x = start
     while x < stop:
@@ -33,6 +36,6 @@ if __name__ == "__main__":
     err = 1e-4
     N0 = 100
 
-    x0 = FindRootZone(f, 5, 6, 0.01)
+    x0 = FindRootZone(f, 5, 6, 0.1) # 查找f的根的区间与步长，改成题干对应的范围
     result,status = NewtonSolve(f, df, x0, err,N0)
     print(f"sqrt(30) = {result:.3f}")

159-7.py

@@ -1,3 +1,4 @@
+# 查找根区间
 def FindRootZone(fx,start,stop,step):
     x = start
     while x < stop:
@@ -32,5 +33,5 @@ def GetNthRoot(a,n):
     return x1
 
 if __name__ == "__main__":
-    re = GetNthRoot(30, 5)
+    re = GetNthRoot(30, 5)   #30的5次方根###########################
     print(re)

160-11.py

@@ -17,11 +17,11 @@ def SecantSolve(fx, x0, x1, err=1e-10, N0=100):
         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

160-12.py

@@ -45,12 +45,12 @@ def MullerSolve(fx,x0,x1,x2,err1,err2,N):
         f2 = f3
         q = h1
 
-
+#精度要求#########################（152页）
 if __name__ == "__main__":
     err1 = 1e-5
     err2 = 1e-5
     N = 100
-
+##把初始函数和初始值改成题干要求的##########################
     x0 = 0.3
     x1 = 0.5
     x2 = 0.4

227-3.py

@@ -19,6 +19,10 @@ def SovleRowMain(A,b):
         b[i], b[row_max_index] = b[row_max_index], b[i]
         p[i], p[row_max_index] = p[row_max_index], p[i]
 
+        print(f"第{i+1}次交换，交换行{i+1}和行{row_max_index+1}，A矩阵为：")
+        for row in A:
+            print(row)
+        print(f"b向量为：{b}\n")
         if abs(A[i][i]) < ks:
             raise ValueError("A矩阵奇异，无法进行高斯消元")
         for j in range(i + 1, n):
@@ -27,7 +31,11 @@ def SovleRowMain(A,b):
             for k in range(i + 1, n):
                 A[j][k] -= m * A[i][k]
             b[j] -= m * b[i]
-
+            print(f"系数为{-1*m}用加号")
+        print("消元后的A矩阵：")
+        for row in A:
+            print(row)
+        print(f"消元后的b向量：{b}\n")
     if abs(A[n - 1][n - 1]) < ks:
         raise ValueError("A矩阵奇异，无法进行高斯消元")
     
@@ -63,7 +71,7 @@ def prettyPrintMatrix(matrix):
 if __name__ == "__main__":
     # A = np.array(A, dtype=float)
     # b = np.array(b, dtype=float)
-      
+       #把矩阵A和b改成题干要求的#####################################
     A = [
         [0, 3, 4],
         [1, -1, 1],

227-4.py

@@ -22,19 +22,31 @@ def SqrtSolve(A,b):
                 for k in range(j):
                     L[getIndexFromDownMatrix(i,j)] -= L[getIndexFromDownMatrix(i,k)]*L[getIndexFromDownMatrix(j,k)]
                 L[getIndexFromDownMatrix(i,j)] /= L[getIndexFromDownMatrix(j,j)]
-
+    # 打印下三角矩阵
+    print("下三角矩阵 L:")
+    for i in range(n):
+        L_row = []
+        for j in range(n):
+            if j <= i:
+                L_row.append(L[getIndexFromDownMatrix(i,j)])
+            else:
+                L_row.append(0)
+        print(L_row)
+    # print(L)
     for i in range(n):
         for k in range(i):
             b[i] -= L[getIndexFromDownMatrix(i,k)]*b[k]
         b[i] /= L[getIndexFromDownMatrix(i,i)]
-
+    # 打印 b 向量
+    print("y 向量:")
+    print(b)
     for i in range(n-1,-1,-1):
         for k in range(i+1,n):
             b[i] -= L[getIndexFromDownMatrix(k,i)]*b[k]
         b[i] /= L[getIndexFromDownMatrix(i,i)]
     return b
 
-
+#把A,b换成题干的数值###########################################
 if __name__ == "__main__":
     # 储存下三角矩阵 a11, a21, a22, a31, a32, a33 ...
     A = [4,2,2,-2,-3,14]
@@ -42,5 +54,5 @@ if __name__ == "__main__":
 
     b = [10,5,4]
 
-    print("x:")
+    # print("x:")
-    print(SqrtSolve(A,b))
+    print("x: \n",SqrtSolve(A,b))

227-7.py

@@ -8,19 +8,21 @@ def ZGsolve(A,b):
             beta[i] = A[i][2] / A[i][1]
         else:
             beta[i] = A[i][2] / (A[i][1] - A[i][0]*beta[i-1])
-
+    print("beta:")
+    print(beta[:-1])
     for i in range(n):
         if i == 0:
             b[i] = b[i] / A[i][1]
         else:
             b[i] = (b[i] - A[i][0]*b[i-1]) / (A[i][1] - A[i][0]*beta[i-1])
-
+    print("y:")
+    print(b)
     for i in range(n-2,-1,-1):
         b[i] = b[i] - beta[i]*b[i+1]
 
     return b
 
-
+#把A,b换成题干的数值###########################################
 if __name__ == "__main__":
     # 储存追赶法A矩阵
     A = [

227-8.py

@@ -19,7 +19,7 @@ def Dot(A,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]]

228-12.py

@@ -51,6 +51,7 @@ def Det(A):
         det += ((-1) ** c) * A[0][c] * Det(sub_matrix)
     return det
 
+#把矩阵换成题干的矩阵#########################
 if __name__ == "__main__":
     A =[
         [1,0,1],

228-13.py

@@ -32,7 +32,7 @@ def SOR(A,b,x,w,err,N):
         
         if count > N:
             return None,count, 0
-        
+#把矩阵改成题干的矩阵，b改成题干结果，err精度要求修改##########################
 if __name__ == "__main__":
     A = [
         [4,-1,0,-1,0,0],
@@ -46,7 +46,7 @@ if __name__ == "__main__":
 
     x = [0,0, 0, 0, 0, 0]
     err = 1e-5
-
+#w换成题干要求的值###########################
     w = 1
     x1,k,sta = SOR(A, b, x, w, err, 100)
     print(f"w = {w}, 解为: {x1}, 迭代次数: {k}, 状态: {'收敛' if sta == 1 else '未收敛'}")

228-15.py

@@ -33,29 +33,51 @@ def Inverse(A):
             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)
+    print("det(A):",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 FixInv(A):
+    n = len(A)
+    if n == 2:
+        return [[A[1][1] / Det(A), -A[0][1] / Det(A)],
+                [-A[1][0] / Det(A), A[0][0] / Det(A)]]
+    B = [[0 for i in range(n)] for j in range(n)]
+    for i in range(n):
+        for j in range(n):
+            t = [row[:j] + row[j+1:] for row in (A[:i] + A[i+1:])]
+            B[j][i] = ((-1) ** (i + j)) * Det(t)
+    det = Det(A)
+    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)
+    # inv_A = Inverse(A)
+    inv_A = FixInv(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]
     ]
-
+#把范数的种类数换成题干的要求，inf是无穷范数#########################
     print(f"矩阵A的条件数为: {Cond(A, float('inf')):.5f}")
-
+#把矩阵换成题干的矩阵#########################
     A = [
         [1,2],
         [3,4]
     ]
+#把范数的种类数换成题干的要求，inf是无穷范数########################
     print(f"矩阵A的条件数为: {Cond(A, float('inf')):.5f}")

228-17.py

@@ -71,14 +71,14 @@ def IterativeMethod(A, b, err, N):
         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}")
+        print(f"第{count+1}次迭代, r{count+1} = {r1}, d{count+1} = {d1}, x{count+2} = {x0}")
         err_now = max(abs(r1[i]) for i in range(len(r1)))
         count += 1
         if err_now < err:
             break
     return x0,count
 
-
+#把矩阵换成题干的矩阵，b换成题干结果，err精度要求修改##########################
 if __name__ == "__main__":
     A = [
         [51,82],
@@ -91,5 +91,5 @@ if __name__ == "__main__":
 
     x = IterativeMethod(A, b, err, N)[0]
     print(f"解为: {x}")
-
+#先用15题的代码求范数与逆矩阵
 

279-3.py

@@ -7,20 +7,20 @@ def ClassicRK(x0,y0,h,xk,fxy):
         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
+        y0 += h*(k1+2*k2+2*k3+k4)/6
         x0 += h
         result.append((x0,y0))
     return result
 
-
+#把下面参数换成题干的#################
 if __name__=="__main__":
-    x0 = 0
+    x0 = 0   #x的左边界换成题干里面的#################
-    y0 = -1
+    y0 = -1  #y的初始值换成题干里面的#################
-    fxy = lambda x,y: x + y
+    fxy = lambda x,y: x + y    #f(x,y)换成题干里面的#################
-    real_fx = lambda x: -x-1
+    real_fx = lambda x: -x-1     #真实函数换成题干里面的#################
-    h = 0.1
+    h = 0.1         #步长换成题干里面的#################
-    xk = 2
+    xk = 2           #x的右边界换成题干里面的#################
     result = ClassicRK(x0, y0, h, xk, fxy)
-    print("x\ty\t\treal_y")
+    print("x\ty\t\treal_y\t\t\t误差")
     for x, y in result:
-        print(f"{x:.2f}\t{y}\t{real_fx(x)}")
+        print(f"{x:.2f}\t{y:.10f}\t\t{real_fx(x):.10f}\t\t\t{abs(y - real_fx(x))}")

68-1.py

@@ -1,18 +1,15 @@
-# 
 
 import math
 
-list_x = [10,11,12,13]
-list_y = [2.3026,2.3979,2.4849,2.5649]
-
 # 定义原函数和其导函数计算结果
 def FxDiff_n(x,n):
     result = 0
     if n == 0:
+        # 下面改成原函数 ####################3########################################
         result = math.log(x)
     else:
+        # 下面改成n阶导数 ##############################################################################
         result = (-1)**(n+1) * math.factorial(n-1) / (x**n)
-
     return result
 
 # 获取与待求x最接近的两个点
@@ -72,8 +69,11 @@ def LagrangeInterpolation(x,list_x,list_y):
         result += temp * list_y[i]
     return result
 
+
+list_x = [10,11,12,13]  #改成题干的数值#############################
+list_y = [2.3026,2.3979,2.4849,2.5649]  #改成题干的数值#############################
 if __name__ == "__main__":
-    print("线性插值 ln11.75 结果为%f, 截断误差%f" % LinearInterpolation(11.75,list_x,list_y))
+    print("线性插值 ln11.75 结果为%f, 截断误差%f" % LinearInterpolation(11.75,list_x,list_y))  #改成题干的数值#############################
-    print("抛物线插值 ln11.75 结果为%f, 截断误差%f" % ParabolaInterpolation(11.75,list_x,list_y))
+    print("抛物线插值 ln11.75 结果为%f, 截断误差%f" % ParabolaInterpolation(11.75,list_x,list_y))  #改成题干的数值#############################
 
     # print(LagrangeInterpolation(11.75,list_x,list_y))

69-3.py

@@ -1,8 +1,5 @@
 import math
 
-list_x = [0.0,0.2,0.4,0.6,0.8]
-list_y = [1.0,1.2214,1.4918,1.8221,2.2255]
-
 # 定义原函数和其导函数计算结果
 def FxDiff_n(x,n):
     return math.exp(x)
@@ -22,15 +19,16 @@ def GetDyList(list_y):
 def NewtonForwardRegression(x,n,list_x):
     h = list_x[1] - list_x[0]
     t = (x - list_x[0]) / h
-    ks = max([abs(FxDiff_n(list_x[i],n+1)) for i in range(n+1)])
+    ks = max([abs(FxDiff_n(list_x[i],n)) for i in range(n)])
     result = 1
-    for i in range(n+1):
+    for i in range(n):
         result *= (t-i)*h/(i+1)
     return  result * ks
 
 # 牛顿前插
 def NewtonForwardInterpolation(x,n,list_x,list_y):
     list_dyk = GetDyList(list_y)
+    print(list_dyk)
     h = list_x[1] - list_x[0]
     t = (x - list_x[0]) / h
     result = list_y[0]
@@ -57,6 +55,7 @@ def GetDQyList(list_x,list_y):
 # 牛顿基本插值
 def NewtonBaseInterpolation(x,n,list_x,list_y):
     list_dqy = GetDQyList(list_x,list_y)
+    print(list_dqy)
     result = list_dqy[0][0]
     mul = 1
     for i in range(1,n):
@@ -65,7 +64,11 @@ def NewtonBaseInterpolation(x,n,list_x,list_y):
     return result
 
 if __name__ == '__main__':
-    print("三点牛顿前插 e^0.12 结果为%f, 截断误差%f" % NewtonForwardInterpolation(0.12, 2, list_x, list_y))
+    
-    print("四点牛顿前插 e^0.12 结果为%f, 截断误差%f" % NewtonForwardInterpolation(0.12, 3, list_x, list_y))
+    list_x = [0.0,0.2,0.4,0.6,0.8]  #改成题干的数值##############################3
-    print("三点牛顿基本插值 e^0.12 结果为%f" % NewtonBaseInterpolation(0.12, 2, list_x, list_y))
+    list_y = [1.0,1.2214,1.4918,1.8221,2.2255]  #改成题干的数值##############################3
-    print("四点牛顿基本插值 e^0.12 结果为%f" % NewtonBaseInterpolation(0.12, 3, list_x, list_y))
+    print("三点牛顿前插 e^0.12 结果为%f, 截断误差%f" % NewtonForwardInterpolation(0.12, 3, list_x, list_y))
+    print("四点牛顿前插 e^0.12 结果为%f, 截断误差%f" % NewtonForwardInterpolation(0.12, 4, list_x, list_y))
+    print("三点牛顿基本插值 e^0.12 结果为%f" % NewtonBaseInterpolation(0.12, 3, list_x, list_y))
+    print("四点牛顿基本插值 e^0.12 结果为%f" % NewtonBaseInterpolation(0.12, 4, list_x, list_y))
+    #70到73行，0.12改成题干的数值，n是点数################################

69-5 - 副本.py

@@ -1,92 +0,0 @@
-x = [0,1,2,3]
-y = [0,0,0,0]
-
-def ZGsolve(A,b):
-    n = len(b)
-    beta = [0]*n
-    for i in range(n):
-        if i == 0:
-            beta[i] = A[i][2] / A[i][1]
-        else:
-            beta[i] = A[i][2] / (A[i][1] - A[i][0]*beta[i-1])
-
-    for i in range(n):
-        if i == 0:
-            b[i] = b[i] / A[i][1]
-        else:
-            b[i] = (b[i] - A[i][0]*b[i-1]) / (A[i][1] - A[i][0]*beta[i-1])
-
-    for i in range(n-2,-1,-1):
-        b[i] = b[i] - beta[i]*b[i+1]
-
-    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)
-    list_mu = [list_h[i]/(list_h[i]+list_h[i+1]) for i in range(len(list_h)-1)]
-    list_lamda = [1-i for i in list_mu]
-    list_g = [6*(list_dqxy[i+1]-list_dqxy[i])/(list_h[i+1]+list_h[i]) for i in range(len(list_h)-1)]
-
-    A = []
-    b = []
-    if boundary_type == 0:  # 自然边界条件
-        pass
-    elif boundary_type == 1:  # 一阶导数边界条件
-        A.append([0,2,1])
-        b.append(6/list_h[0]*(list_dqxy[1]-a1))
-        for i in range(len(list_g)):
-            A.append([list_mu[i],2,list_lamda[i]])
-            b.append(list_g[i])
-        A.append([1,2,0])
-        b.append(6/list_h[-1]*(a2-list_dqxy[-1]))
-        M = ZGsolve(A,b)
-
-    elif boundary_type == 2:  # 二阶导数边界条件
-        A.append([0,2,list_lamda[0]])
-        b.append(list_g[0]-list_mu[0]*a1)
-        for i in range(1,len(list_g)-1):
-            A.append([list_mu[i],2,list_lamda[i]])
-            b.append(list_g[i])
-        A.append([list_mu[-1],2,0])
-        b.append(list_g[-1]-list_lamda[-1]*a2)
-        M = ZGsolve(A,b)
-        M = [a1] + M + [a2]
-
-    return M,list_h
-
-
-if __name__ == "__main__":
-    M,list_h = CubicSplineInterpolation(x,y,2,1,0)
-
-    print("list_h:",list_h)
-    print("M:",M)
-    for i in range(3):
-        k1 = (M[i+1]-M[i])/6/list_h[i]
-        k2 = (M[i]*x[i+1]-M[i+1]*x[i])/2/list_h[i]
-        k3 = (3*M[i+1]*x[i]**2-3*M[i]*x[i+1]**2-6*y[i]+M[i]*list_h[i]**2+6*y[i+1]-M[i+1]*list_h[i]**2)/6/list_h[i]
-        k4 = (M[i]*x[i+1]**3-M[i+1]*x[i]**3+6*y[i]*x[i+1]-M[i]*list_h[i]**2*x[i+1]-6*y[i+1]*x[i]+M[i+1]*list_h[i]**2*x[i])/6/list_h[i]
-        # print(k1,k2,k3,k4)
-        print("S(x)=%.10f*x^3+%.10f*x^2+%.10f*x+%.10f"%(k1,k2,k3,k4))
-
-
-
-
-
-    
-
-
-

69-5.py

@@ -1,6 +1,5 @@
-x = [0,1,2,3]
-y = [0,0,0,0]
 
+# 追赶法
 def ZGsolve(A,b):
     n = len(b)
     beta = [0]*n
@@ -21,30 +20,50 @@ 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)
+    print("h:", list_h)
     list_dqxy = GetDQList(list_x, list_y)
+    print("f[xi,xi+1]:", list_dqxy)
     list_mu = [list_h[i]/(list_h[i]+list_h[i+1]) for i in range(len(list_h)-1)]
+    print("miu:", list_mu)
     list_lamda = [1-i for i in list_mu]
+    print("lambda:", list_lamda)
     list_g = [6*(list_dqxy[i+1]-list_dqxy[i])/(list_h[i+1]+list_h[i]) for i in range(len(list_h)-1)]
     
     A = []
     b = []
+    M = []
+    copy_b = []
     if boundary_type == 0:  # 自然边界条件
-        pass
+        a1 = 0
+        a2 = 0
+        A.append([0,2,list_lamda[0]])
+        b.append(list_g[0]-list_mu[0]*a1)
+        for i in range(1,len(list_g)-1):
+            A.append([list_mu[i],2,list_lamda[i]])
+            b.append(list_g[i])
+        A.append([list_mu[-1],2,0])
+        b.append(list_g[-1]-list_lamda[-1]*a2)
+        copy_b = b.copy()
+        print("g1~gn-1:", list_g)
+        M = ZGsolve(A,b)
+        M = [a1] + M + [a2]
     elif boundary_type == 1:  # 一阶导数边界条件
         A.append([0,2,1])
         b.append(6/list_h[0]*(list_dqxy[1]-a1))
@@ -53,6 +72,9 @@ def CubicSplineInterpolation(list_x,list_y,boundary_type,a1,a2):
             b.append(list_g[i])
         A.append([1,2,0])
         b.append(6/list_h[-1]*(a2-list_dqxy[-1]))
+        copy_b = b.copy()
+        print("g0~gn:", copy_b)
+        
         M = ZGsolve(A,b)
 
     elif boundary_type == 2:  # 二阶导数边界条件
@@ -63,11 +85,16 @@ def CubicSplineInterpolation(list_x,list_y,boundary_type,a1,a2):
             b.append(list_g[i])
         A.append([list_mu[-1],2,0])
         b.append(list_g[-1]-list_lamda[-1]*a2)
+        copy_b = b.copy()
+        print("g1~gn-1:", list_g)
         M = ZGsolve(A,b)
         M = [a1] + M + [a2]
-
+    print("A:", A)
+    print("b:", copy_b)
+    print("M:", M)
     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]
@@ -76,8 +103,10 @@ def PrintResult(M,list_h,list_x,list_y):
         k4 = (M[i]*list_x[i+1]**3-M[i+1]*list_x[i]**3+6*list_y[i]*list_x[i+1]-M[i]*list_h[i]**2*list_x[i+1]-6*list_y[i+1]*list_x[i]+M[i+1]*list_h[i]**2*list_x[i])/6/list_h[i]
         print("S(x)=%.6f*x^3+%.6f*x^2+%6f*x+%6f"%(k1,k2,k3,k4),"x=[%.6f,%.6f]"%(list_x[i],list_x[i+1]))
 
-
+#将x,y换成题干的数值###########################################
 if __name__ == "__main__":
+    x = [0,1,2,3]
+    y = [0,0,0,0]
     M,list_h = CubicSplineInterpolation(x,y,2,1,0)
     print("二阶导数边界条件:")
     PrintResult(M,list_h,x,y)

89-1.py

@@ -1,6 +1,3 @@
-x = [2,4,6,8]
-y = [2,11,28,40]
-
 
 # 列主元高斯消元法
 def SovleRowMain(A,b):
@@ -43,9 +40,17 @@ def SovleRowMain(A,b):
     return b
 
 
-
+# 最小二乘法拟合
 def LeastSquares(list_x,list_y,n):
     m = len(list_x)
+    for i in range(2*n):
+        tl = [list_x[j]**(i+1) for j in range(m)]
+        print(f"x^{i+1} :",tl)
+        print("sum:", sum(tl))
+    for i in range(n+1):
+        tl = [list_y[j]*list_x[j]**i for j in range(m)]
+        print(f"y*x^{i}: ",tl)
+        print("sum:", sum(tl))
     x_n = []
     b = []
     for i in range(2*n+1):
@@ -61,23 +66,33 @@ def LeastSquares(list_x,list_y,n):
         for j in range(n+1):
             tmp.append(x_n[i+j])
         A.append(tmp)
-    return SovleRowMain(A,b)
+    print("A:", A)
+    print("b:", b)
+    result = SovleRowMain(A, b)
+    print("result:", result)
+    return result
 
+# 计算多项式在给定x值上的值
 def CalculateY(list_x, coeff):
     re = []
     for i in range(len(list_x)):
         re.append(0)
         for j in range(len(coeff)):
             re[i] += coeff[j]*list_x[i]**j
+    print("y预测:", re)
     return re
 
+# 计算均方根误差
 def MeanSquareErr(list_y,list_y_approx):
     m = len(list_y)
     err = 0
+    print("y差值: ")
     for i in range(m):
+        print(f"y-y测={list_y[i] - list_y_approx[i]}, (y-y测)^2={(list_y[i] - list_y_approx[i])**2}")
         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 +101,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]) + "+"
@@ -95,22 +111,27 @@ def PrintEquation(coeff):
     str_ = str_[0:len(str_)-1]
     print(str_)
 
+#把x和y换成题干的数值###################
+x = [2,4,6,8]
+y = [2,11,28,40]
 if __name__ == "__main__":
     print("一次拟合")
     coeff = LeastSquares(x,y,1)
     PrintEquation(coeff)
     y_approx = CalculateY(x, coeff)
-    print("MeanSquareErr:")
+    # print()
-    print(MeanSquareErr(y,y_approx))
+    print("均方根误差:",MeanSquareErr(y,y_approx))
-    print("MaxErr:")
+    # print()
-    print(MaxErr(y,y_approx))
+    print("最大误差:",MaxErr(y,y_approx))
+
+    print()
 
     print("二次拟合")
     coeff = LeastSquares(x,y,2)
     PrintEquation(coeff)
     y_approx = CalculateY(x, coeff)
-    print("MeanSquareErr:")
+    # print()
-    print(MeanSquareErr(y,y_approx))
+    print("均方根误差:",MeanSquareErr(y,y_approx))
-    print("MaxErr:")
+    # print()
-    print(MaxErr(y,y_approx))
+    print("最大误差:",MaxErr(y,y_approx))
     

89-2 - 副本.py

@@ -1,76 +0,0 @@
-import math
-
-x = [1,2,4,8,16,32,64]
-y = [4.22,4.02,3.85,3.59,3.44,3.02,2.59]
-
-
-# 列主元高斯消元法
-def SovleRowMain(A,b):
-    ks = 0.00000001
-    n = len(A)
-    if len(A[0]) != n:
-        raise ValueError("A要为方阵")
-    if len(b) != n:
-        raise ValueError("b与A的行数不匹配")
-    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:
-            raise ValueError("A矩阵奇异，无法进行高斯消元")
-        for j in range(i + 1, n):
-            m = A[j][i] / A[i][i]
-            A[j][i] = m
-            for k in range(i + 1, n):
-                A[j][k] -= m * A[i][k]
-            b[j] -= m * b[i]
-
-    if abs(A[n - 1][n - 1]) < ks:
-        raise ValueError("A矩阵奇异，无法进行高斯消元")
-    
-    # 回代求解
-    b[n - 1] /= A[n - 1][n - 1]
-    for i in range(n - 2, -1, -1):
-        for j in range(i + 1, n):
-            b[i] -= A[i][j] * b[j]
-        b[i] /= A[i][i]
-    return b
-
-
-
-def LeastSquares(list_x,list_y,n):
-    m = len(list_x)
-    x_n = []
-    b = []
-    for i in range(2*n+1):
-        x_n.append(0)
-        b.append(0)
-        for j in range(m):
-            x_n[i]+=(list_x[j]**i)
-            b[i]+=(list_y[j]*list_x[j]**i)
-    b = b[:n+1]
-    A = []
-    for i in range(n+1):
-        tmp = []
-        for j in range(n+1):
-            tmp.append(x_n[i+j])
-        A.append(tmp)
-    return SovleRowMain(A,b)
-
-if __name__ == "__main__":
-    # 取对数 ln(W) = ln(C)+lamda*ln(t)
-    ln_W = [math.log(i) for i in y]
-    ln_t = [math.log(i) for i in x]
-    coeff = LeastSquares(ln_t,ln_W,1)
-    C = math.exp(coeff[0])
-    lamda = coeff[1]
-    print("C:", C)
-    print("lamda:", lamda)

89-2.py

@@ -1,9 +1,5 @@
 import math
 
-x = [1,2,4,8,16,32,64]
-y = [4.22,4.02,3.85,3.59,3.44,3.02,2.59]
-
-
 # 列主元高斯消元法
 def SovleRowMain(A,b):
     ks = 0.00000001
@@ -45,7 +41,7 @@ def SovleRowMain(A,b):
     return b
 
 
-
+# 最小二乘法拟合
 def LeastSquares(list_x,list_y,n):
     m = len(list_x)
     x_n = []
@@ -63,8 +59,15 @@ def LeastSquares(list_x,list_y,n):
         for j in range(n+1):
             tmp.append(x_n[i+j])
         A.append(tmp)
-    return SovleRowMain(A,b)
+    print("A:", A)
+    print("b:", b)
+    result = SovleRowMain(A, b)
+    print("result:", result)
+    return result
 
+#把x和y的值改为实际数据#######################
+x = [1,2,4,8,16,32,64]
+y = [4.22,4.02,3.85,3.59,3.44,3.02,2.59]
 if __name__ == "__main__":
     # 取对数 ln(W) = ln(C)+lamda*ln(t)
     ln_W = [math.log(i) for i in y]

89-3 - 副本.py

@@ -1,69 +0,0 @@
-
-x = [19,25,31,38,44]
-y = [19.0,32.3,49.0,73.3,97.8]
-
-# 列主元高斯消元法
-def SovleRowMain(A,b):
-    ks = 0.00000001
-    n = len(A)
-    if len(A[0]) != n:
-        raise ValueError("A要为方阵")
-    if len(b) != n:
-        raise ValueError("b与A的行数不匹配")
-    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:
-            raise ValueError("A矩阵奇异，无法进行高斯消元")
-        for j in range(i + 1, n):
-            m = A[j][i] / A[i][i]
-            A[j][i] = m
-            for k in range(i + 1, n):
-                A[j][k] -= m * A[i][k]
-            b[j] -= m * b[i]
-
-    if abs(A[n - 1][n - 1]) < ks:
-        raise ValueError("A矩阵奇异，无法进行高斯消元")
-    
-    # 回代求解
-    b[n - 1] /= A[n - 1][n - 1]
-    for i in range(n - 2, -1, -1):
-        for j in range(i + 1, n):
-            b[i] -= A[i][j] * b[j]
-        b[i] /= A[i][i]
-    return b
-
-
-
-def LeastSquares(list_x,list_y,n):
-    m = len(list_x)
-    x_n = []
-    b = []
-    for i in range(2*n+1):
-        x_n.append(0)
-        b.append(0)
-        for j in range(m):
-            x_n[i]+=(list_x[j]**i)
-            b[i]+=(list_y[j]*list_x[j]**i)
-    b = b[:n+1]
-    A = []
-    for i in range(n+1):
-        tmp = []
-        for j in range(n+1):
-            tmp.append(x_n[i+j])
-        A.append(tmp)
-    return SovleRowMain(A,b)
-
-if __name__ == "__main__":
-    x_square = [i**2 for i in x]
-    coeff = LeastSquares(x_square, y, 1)
-    print("拟合方程: y = %.6f + %.6f*x^2" % (coeff[0], coeff[1]))

89-3.py

@@ -1,7 +1,4 @@
 
-x = [19,25,31,38,44]
-y = [19.0,32.3,49.0,73.3,97.8]
-
 # 列主元高斯消元法
 def SovleRowMain(A,b):
     ks = 0.00000001
@@ -43,7 +40,7 @@ def SovleRowMain(A,b):
     return b
 
 
-
+# 最小二乘法拟合
 def LeastSquares(list_x,list_y,n):
     m = len(list_x)
     x_n = []
@@ -61,8 +58,15 @@ def LeastSquares(list_x,list_y,n):
         for j in range(n+1):
             tmp.append(x_n[i+j])
         A.append(tmp)
-    return SovleRowMain(A,b)
+    print("A:", A)
+    print("b:", b)
+    result = SovleRowMain(A, b)
+    print("result:", result)
+    return result
 
+#把x和y换成题干的数值###################
+x = [19,25,31,38,44]
+y = [19.0,32.3,49.0,73.3,97.8]
 if __name__ == "__main__":
     x_square = [i**2 for i in x]
     coeff = LeastSquares(x_square, y, 1)

按方法整理/常微分方程-经典龙格-库塔格式.py

@@ -0,0 +1,29 @@
+import math
+
+#经典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+2*k2+2*k3+k4)/6
+        x0 += h
+        result.append((x0,y0))
+    return result
+
+
+if __name__=="__main__":
+    ##########################################################################
+    x0 = 0   # x的初始值（左边界）换成题干里面的#################
+    y0 = -1  # y的初始值换成题干里面的#################
+    fxy = lambda x,y: x+y   #f(x,y)换成题干里面的#################
+    real_fx = lambda x: -x-1     #真实函数换成题干里面的#################
+    h = 0.1         #步长换成题干里面的#################
+    xk = 2           #x的右边界换成题干里面的#################
+    result = ClassicRK(x0, y0, h, xk, fxy)
+    print("x\ty\t\treal_y\t\t\t误差")
+    for x, y in result:
+        print(f"{x:.2f}\t{y:.10f}\t\t{real_fx(x):.10f}\t\t\t{abs(y - real_fx(x)):.5f}")

按方法整理/常微分方程-阿达姆斯方法.py

@@ -0,0 +1,54 @@
+import math
+
+
+def AdamusExplicitly(k,first_ys,x0,y0,h,xk,fxy,fx_real):
+    betas = [
+        [1],
+        [3/2, -1/2],
+        [23/12, -16/12, 5/12],
+        [55/24, -59/24, 37/24, -9/24],
+        [1901/720, -2774/720, 2626/720, -1274/720, 251/720],
+        [4277/1440, -7923/1600, 9982/1440, -7298/1440, 2877/1440, -476/1440]
+    ]
+    # Bks = [1/2,5/12,3/8,251/720,95/288,10987/60480]
+    result = []
+    if k<0 or k >= len(betas):
+        print("k超出范围")
+        return None
+    if len(first_ys) != k + 1:
+        print("first_ys(前几个y的值)与k长度不匹配")
+        return None
+    fxys = [fxy(x0+i*h, first_ys[i]) for i in range(k + 1)]
+    for i in range(k + 1):
+        result.append((x0 + i * h, first_ys[i]))
+    x1 = x0 + h*k
+    y0 = first_ys[k]
+    while x1 < xk:
+        x1 += h
+        delta_y = h * sum(betas[k][i] * fxys[k-i] for i in range(k + 1))
+        y1 = y0 + delta_y
+        print(f"x={x1}, y={y1}, delta_y={delta_y}, real_y={fx_real(x1)}, abs(real_y-y)={abs(fx_real(x1) - y1)}")
+        result.append((x1, y1))
+        fxys.append(fxy(x1, y1))
+        fxys.pop(0)
+        y0 = y1
+    return result
+
+
+
+if __name__ == "__main__":
+    # 定义初始值和参数#####################################################################################
+    k = 3                    # 阿达姆斯k+1步显式方法，精度为k+1阶，最常用k=3，其他阶数我没试过   P253
+    first_ys = [1, 0.904837418036, 0.818730753078, 0.7408182206817]  # 前几个y的值，可用龙格-库塔计算或者知道精确解自己算出来
+    x0 = 0.0                 # 初始x值
+    y0 = 1.0                 # 初始y值
+    h = 0.1                  # 步长
+    xk = 1.0                 # 最终x值
+    fxy = lambda x, y: -y                   # 定义f(x,y)函数，dx/dy = f(x,y)，导函数
+    fx_real = lambda x: math.exp(-x)        # 实际解函数，用于验证结果，如果不知道或者不用算误差，可以直接写个 lambda x: 0
+    # 调用阿达姆斯显式方法
+    result = AdamusExplicitly(k, first_ys, x0, y0, h, xk, fxy, fx_real)
+    print("计算结果看到有几.99999或者几.00000就自己四舍五入一下，有可能会多算一点，自己比较一下")
+    if result:
+        for xy in result:
+            print(f"(x, y): {xy}")

按方法整理/常微分方程-阿达姆斯隐式方法.py

@@ -0,0 +1,61 @@
+import math
+
+
+def AdamusImplicitly(k,first_ys,x0,y0,h,xk,fxy,fx_real):
+    betas = [
+        [1],
+        [1/2, 1/2],
+        [5/12, 8/12, -1/12],
+        [9/24, 19/24, -5/24, 1/24],
+        [251/720, 646/720, -261/720, 106/720, -19/720],
+        [475/1440, 1427/1440, -798/1440, 482/1440, -173/1440, 27/1440]
+    ]
+    
+    result = []
+    if k<0 or k >= len(betas):
+        print("k超出范围")
+        return None
+    if len(first_ys) != k:
+        print("first_ys(前几个y的值)与k长度不匹配")
+        return None
+    fxys = [fxy(x0+i*h, first_ys[i]) for i in range(k)]
+    for i in range(k):
+        result.append((x0 + i * h, first_ys[i]))
+    x1 = x0 + h*k
+    y0 = first_ys[k-1]
+    while x1 < xk:
+        y1_t = y0
+        y1 = y0 + 1
+        for j in range(1000):
+            y1 = y0 + h * sum(betas[k][i+1] * fxys[k-i-1] for i in range(k))+ h * betas[k][0] * fxy(x1, y1_t)
+            if abs(y1 - y1_t) < 1e-14:
+                break
+            y1_t = y1
+
+        # y1 = y0 + delta_y
+        print(f"x={x1}, y={y1}, real_y={fx_real(x1)}, abs(real_y-y)={abs(fx_real(x1) - y1)}")
+        result.append((x1, y1))
+        fxys.append(fxy(x1, y1))
+        fxys.pop(0)
+        y0 = y1
+        x1 += h
+    return result
+
+
+
+if __name__ == "__main__":
+    # 定义初始值和参数#####################################################################################
+    k = 3                    # 阿达姆斯k步隐式方法，精度为k+1阶，最常用k=3，其他阶数我没试过   P255
+    first_ys = [1, 0.904837418036, 0.818730753078]  # 前几个y的值，可用龙格-库塔计算或者知道精确解自己算出来
+    x0 = 0.0                 # 初始x值
+    y0 = 1.0                 # 初始y值
+    h = 0.1                  # 步长
+    xk = 1.0                 # 最终x值
+    fxy = lambda x, y: -y                   # 定义f(x,y)函数，dx/dy = f(x,y)，导函数
+    fx_real = lambda x: math.exp(-x)        # 实际解函数，用于验证结果，如果不知道或者不用算误差，可以直接写个 lambda x: 0
+    # 调用阿达姆斯隐式方法
+    result = AdamusImplicitly(k, first_ys, x0, y0, h, xk, fxy, fx_real)
+    print("计算结果看到有几.99999或者几.00000就自己四舍五入一下，有可能会多算一点，自己比较一下")
+    if result:
+        for xy in result:
+            print(f"(x, y): {xy}")

按方法整理/常微分方程-阿达姆斯预测-校正方法.py

@@ -0,0 +1,65 @@
+import math
+
+
+def AdamusFix(k,first_ys,x0,y0,h,xk,fxy,fx_real):
+    betas_1 = [
+        [1],
+        [3/2, -1/2],
+        [23/12, -16/12, 5/12],
+        [55/24, -59/24, 37/24, -9/24],
+        [1901/720, -2774/720, 2626/720, -1274/720, 251/720],
+        [4277/1440, -7923/1600, 9982/1440, -7298/1440, 2877/1440, -476/1440]
+    ]
+    # Bks = [1/2,5/12,3/8,251/720,95/288,10987/60480]
+    betas_2 = [
+        [1],
+        [1/2, 1/2],
+        [5/12, 8/12, -1/12],
+        [9/24, 19/24, -5/24, 1/24],
+        [251/720, 646/720, -261/720, 106/720, -19/720],
+        [475/1440, 1427/1440, -798/1440, 482/1440, -173/1440, 27/1440]
+    ]
+    result = []
+    if k<0 or k >= len(betas_1):
+        print("k超出范围")
+        return None
+    if len(first_ys) != k + 1:
+        print("first_ys(前几个y的值)与k长度不匹配")
+        return None
+    fxys = [fxy(x0+i*h, first_ys[i]) for i in range(k + 1)]
+    for i in range(k + 1):
+        result.append((x0 + i * h, first_ys[i]))
+    x1 = x0 + h*k
+    y0 = first_ys[k]
+    while x1 < xk:
+        x1 += h
+        y1_guess = y0 + h * sum(betas_1[k][i] * fxys[k-i] for i in range(k + 1))
+        fxy_guess = fxy(x1, y1_guess)
+        fxys.append(fxy_guess)
+        fxys.pop(0)
+        delta_y = h * sum(betas_2[k][i] * fxys[k-i] for i in range(k + 1))
+        y1 = y0 + delta_y
+        print(f"x={x1}, y={y1}, y_guess={y1_guess}, real_y={fx_real(x1)}, abs(real_y-y)={abs(fx_real(x1) - y1)}")
+        result.append((x1, y1))
+        fxys[k] = fxy(x1, y1)
+        y0 = y1
+    return result
+
+
+
+if __name__ == "__main__":
+    # 定义初始值和参数#####################################################################################
+    k = 3                    # 阿达姆斯k+1步校正方法，精度为k+1阶，最常用k=3，其他阶数我没试过   P253
+    first_ys = [1, 1.0954, 1.1832, 1.2649]  # 前几个y的值，可用龙格-库塔计算或者知道精确解自己算出来
+    x0 = 0.0                 # 初始x值
+    y0 = 1.0                 # 初始y值
+    h = 0.1                  # 步长
+    xk = 1.0                 # 最终x值
+    fxy = lambda x, y: y - 2*x/y          # 定义f(x,y)函数，dx/dy = f(x,y)，导函数
+    fx_real = lambda x: 0        # 实际解函数，用于验证结果，如果不知道或者不用算误差，可以直接写个 lambda x: 0
+    # 调用阿达姆斯显式方法
+    result = AdamusFix(k, first_ys, x0, y0, h, xk, fxy, fx_real)
+    print("计算结果看到有几.99999或者几.00000就自己四舍五入一下，有可能会多算一点，自己比较一下")
+    if result:
+        for xy in result:
+            print(f"(x, y): {xy}")

按方法整理/拟合-最小二乘法.py

@@ -1,15 +1,15 @@
-x = [2,4,6,8]
+import math
-y = [2,11,28,40]
-
 
 # 列主元高斯消元法
 def SovleRowMain(A,b):
     ks = 0.00000001
     n = len(A)
     if len(A[0]) != n:
-        raise ValueError("A要为方阵")
+        print("A要为方阵")
+        return None,None,None,None
     if len(b) != n:
-        raise ValueError("b与A的行数不匹配")
+        print("b与A的行数不匹配")
+        return None,None,None,None
     p = list(range(n))
     for i in range(n):
         row_max = abs(A[i][i])
@@ -23,7 +23,8 @@ def SovleRowMain(A,b):
         p[i], p[row_max_index] = p[row_max_index], p[i]
 
         if abs(A[i][i]) < ks:
-            raise ValueError("A矩阵奇异，无法进行高斯消元")
+            print("A矩阵奇异，无法进行高斯消元")
+            return None,None,None,None
         for j in range(i + 1, n):
             m = A[j][i] / A[i][i]
             A[j][i] = m
@@ -32,7 +33,8 @@ def SovleRowMain(A,b):
             b[j] -= m * b[i]
 
     if abs(A[n - 1][n - 1]) < ks:
-        raise ValueError("A矩阵奇异，无法进行高斯消元")
+        print("A矩阵奇异，无法进行高斯消元")
+        return None,None,None,None
     
     # 回代求解
     b[n - 1] /= A[n - 1][n - 1]
@@ -43,9 +45,17 @@ def SovleRowMain(A,b):
     return b
 
 
-
+# 最小二乘法拟合
 def LeastSquares(list_x,list_y,n):
     m = len(list_x)
+    for i in range(2*n):
+        tl = [list_x[j]**(i+1) for j in range(m)]
+        print(f"x^{i+1} :",tl)
+        print("sum:", sum(tl))
+    for i in range(n+1):
+        tl = [list_y[j]*list_x[j]**i for j in range(m)]
+        print(f"y*x^{i}: ",tl)
+        print("sum:", sum(tl))
     x_n = []
     b = []
     for i in range(2*n+1):
@@ -61,23 +71,33 @@ def LeastSquares(list_x,list_y,n):
         for j in range(n+1):
             tmp.append(x_n[i+j])
         A.append(tmp)
-    return SovleRowMain(A,b)
+    print("A:", A)
+    print("b:", b)
+    result = SovleRowMain(A, b)
+    print("result:", result)
+    return result
 
+# 计算多项式在给定x值上的值
 def CalculateY(list_x, coeff):
     re = []
     for i in range(len(list_x)):
         re.append(0)
         for j in range(len(coeff)):
             re[i] += coeff[j]*list_x[i]**j
+    print("y预测:", re)
     return re
 
+# 计算均方根误差
 def MeanSquareErr(list_y,list_y_approx):
     m = len(list_y)
     err = 0
+    print("y差值: ")
     for i in range(m):
+        print(f"y-y测={list_y[i] - list_y_approx[i]}, (y-y测)^2={(list_y[i] - list_y_approx[i])**2}")
         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 +106,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]) + "+"
@@ -95,22 +116,28 @@ def PrintEquation(coeff):
     str_ = str_[0:len(str_)-1]
     print(str_)
 
+
 if __name__ == "__main__":
+    ##########################################################################
+    # 下面的数值根据题意改成适合的，比如要取对数就先取对数，下面的x和y是直接用于拟合的
+    list_x = [2,4,6,8]        # 已给出的x数值，与y数值对应
+    list_y = [2,11,28,40]     # 已给出的y数值，与x数值对应
+
+    # 记得改下面最后一个参数，为拟合阶数
     print("一次拟合")
-    coeff = LeastSquares(x,y,1)
+    coeff = LeastSquares(list_x,list_y,1)
     PrintEquation(coeff)
-    y_approx = CalculateY(x, coeff)
+    y_approx = CalculateY(list_x, coeff)
-    print("MeanSquareErr:")
+    print("均方根误差:",MeanSquareErr(list_y,y_approx))
-    print(MeanSquareErr(y,y_approx))
+    print("最大误差:",MaxErr(list_y,y_approx))
-    print("MaxErr:")
+
-    print(MaxErr(y,y_approx))
+    print()
     
     print("二次拟合")
-    coeff = LeastSquares(x,y,2)
+    coeff = LeastSquares(list_x,list_y,2)
     PrintEquation(coeff)
-    y_approx = CalculateY(x, coeff)
+    y_approx = CalculateY(list_x, coeff)
-    print("MeanSquareErr:")
+    print("均方根误差:",MeanSquareErr(list_y,y_approx))
-    print(MeanSquareErr(y,y_approx))
+    print("最大误差:",MaxErr(list_y,y_approx))
-    print("MaxErr:")
+
-    print(MaxErr(y,y_approx))
     

按方法整理/插值-三次样条.py

@@ -0,0 +1,122 @@
+import math
+
+# 追赶法
+def ZGsolve(A,b):
+    n = len(b)
+    beta = [0]*n
+    for i in range(n):
+        if i == 0:
+            beta[i] = A[i][2] / A[i][1]
+        else:
+            beta[i] = A[i][2] / (A[i][1] - A[i][0]*beta[i-1])
+
+    for i in range(n):
+        if i == 0:
+            b[i] = b[i] / A[i][1]
+        else:
+            b[i] = (b[i] - A[i][0]*b[i-1]) / (A[i][1] - A[i][0]*beta[i-1])
+
+    for i in range(n-2,-1,-1):
+        b[i] = b[i] - beta[i]*b[i+1]
+
+    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)
+    print("h:", list_h)
+    list_dqxy = GetDQList(list_x, list_y)
+    print("f[xi,xi+1]:", list_dqxy)
+    list_mu = [list_h[i]/(list_h[i]+list_h[i+1]) for i in range(len(list_h)-1)]
+    print("miu:", list_mu)
+    list_lamda = [1-i for i in list_mu]
+    print("lambda:", list_lamda)
+    list_g = [6*(list_dqxy[i+1]-list_dqxy[i])/(list_h[i+1]+list_h[i]) for i in range(len(list_h)-1)]
+    A = []
+    b = []
+    M = []
+    copy_b = []
+    if boundary_type == 0:  # 自然边界条件
+        a1 = 0
+        a2 = 0
+        A.append([0,2,list_lamda[0]])
+        b.append(list_g[0]-list_mu[0]*a1)
+        for i in range(1,len(list_g)-1):
+            A.append([list_mu[i],2,list_lamda[i]])
+            b.append(list_g[i])
+        A.append([list_mu[-1],2,0])
+        b.append(list_g[-1]-list_lamda[-1]*a2)
+        copy_b = b.copy()
+        print("g1~gn-1:", list_g)
+        M = ZGsolve(A,b)
+        M = [a1] + M + [a2]
+    elif boundary_type == 1:  # 一阶导数边界条件
+        A.append([0,2,1])
+        b.append(6/list_h[0]*(list_dqxy[1]-a1))
+        for i in range(len(list_g)):
+            A.append([list_mu[i],2,list_lamda[i]])
+            b.append(list_g[i])
+        A.append([1,2,0])
+        b.append(6/list_h[-1]*(a2-list_dqxy[-1]))
+        copy_b = b.copy()
+        print("g0~gn:", copy_b)
+        
+        M = ZGsolve(A,b)
+
+    elif boundary_type == 2:  # 二阶导数边界条件
+        A.append([0,2,list_lamda[0]])
+        b.append(list_g[0]-list_mu[0]*a1)
+        for i in range(1,len(list_g)-1):
+            A.append([list_mu[i],2,list_lamda[i]])
+            b.append(list_g[i])
+        A.append([list_mu[-1],2,0])
+        b.append(list_g[-1]-list_lamda[-1]*a2)
+        copy_b = b.copy()
+        print("g1~gn-1:", list_g)
+        M = ZGsolve(A,b)
+        M = [a1] + M + [a2]
+    print("A:", A)
+    print("b:", copy_b)
+    print("M:", M)
+    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]
+        k2 = (M[i]*list_x[i+1]-M[i+1]*list_x[i])/2/list_h[i]
+        k3 = (3*M[i+1]*list_x[i]**2-3*M[i]*list_x[i+1]**2-6*list_y[i]+M[i]*list_h[i]**2+6*list_y[i+1]-M[i+1]*list_h[i]**2)/6/list_h[i]
+        k4 = (M[i]*list_x[i+1]**3-M[i+1]*list_x[i]**3+6*list_y[i]*list_x[i+1]-M[i]*list_h[i]**2*list_x[i+1]-6*list_y[i+1]*list_x[i]+M[i+1]*list_h[i]**2*list_x[i])/6/list_h[i]
+        print("S(x)=%.6f*x^3+%.6f*x^2+%6f*x+%6f"%(k1,k2,k3,k4),"x=[%.6f,%.6f]"%(list_x[i],list_x[i+1]))
+
+
+if __name__ == "__main__":
+    ##############################################################################################################
+    list_x = [0,1,2,3]               # 已给出的x数值，与y数值对应
+    list_y = [0,0,0,0]               # 已给出的y数值，与x数值对应
+    # 记得改后三个参数，分别为边界条件类型(0=自然边界条件，1=一阶导数边界条件，2=二阶导数边界条件)，边界条件a1和a2
+    M,list_h = CubicSplineInterpolation(list_x,list_y,2,1,0)
+    print("二阶导数边界条件:")
+    PrintResult(M,list_h,list_x,list_y)
+    
+    M,list_h = CubicSplineInterpolation(list_x,list_y,1,1,0)
+    print("一阶导数边界条件:")
+    PrintResult(M,list_h,list_x,list_y)
+
+
+
+

按方法整理/插值-牛顿前插-牛顿基本插值.py

@@ -0,0 +1,93 @@
+import math
+
+# 求差分表
+def GetDyList(list_y):
+    result = []
+    result.append(list_y)
+    for i in range(len(list_y)-1,0,-1):
+        tmp = []
+        for j in range(i):
+            tmp.append(result[len(list_y)-1-i][j+1] - result[len(list_y)-1-i][j])
+        result.append(tmp)
+    return result
+
+# 牛顿前插余项计算
+def NewtonForwardRegression(x,n,list_x,FxDiff_n):
+    h = list_x[1] - list_x[0]
+    t = (x - list_x[0]) / h
+    ks = max([abs(FxDiff_n(list_x[i],n)) for i in range(n)])
+    result = 1
+    for i in range(n):
+        result *= (t-i)*h/(i+1)
+    return  result * ks
+
+# 牛顿前插
+def NewtonForwardInterpolation(x,n,list_x,list_y,FxDiff_n):
+    list_dyk = GetDyList(list_y)
+    print("\n差分表")
+    print("--------------------------------------------------")
+    for i in range(len(list_dyk)):
+        print(list_dyk[i])
+    print("--------------------------------------------------")
+    h = list_x[1] - list_x[0]
+    t = (x - list_x[0]) / h
+    result = list_y[0]
+    mul = 1
+    result = list_y[0]
+    for i in range(1, n):
+        mul *= ((t-i+1)/i)
+        result += list_dyk[i][0] * mul
+    r = abs(NewtonForwardRegression(x,n,list_x,FxDiff_n))
+    return (result,r)
+
+# 求差商表
+def GetDQyList(list_x,list_y):
+    result = []
+    result.append(list_y)
+    for i in range(len(list_y)-1,0,-1):
+        tmp = []
+        for j in range(i):
+            tmp.append((result[len(list_y)-1-i][j+1] - result[len(list_y)-1-i][j])/(list_x[j+len(list_y)-i] - list_x[j]))
+        result.append(tmp)
+    return result
+
+
+# 牛顿基本插值
+def NewtonBaseInterpolation(x,n,list_x,list_y):
+    list_dqy = GetDQyList(list_x,list_y)
+    print("\n差商表")
+    print("--------------------------------------------------")
+    for i in range(len(list_dqy)):
+        print(list_dqy[i])
+    print("--------------------------------------------------")
+    result = list_dqy[0][0]
+    mul = 1
+    for i in range(1,n):
+        mul *= (x - list_x[i-1])
+        result += list_dqy[i][0] * mul
+    return result
+
+
+
+# 定义原函数和其导函数计算结果，用于计算牛顿前插的截断误差，如果不需要则不用管
+def FxDiff_n1(x,n):
+    result = 0
+    if n == 0:
+        # 下面改成原函数 ############################################################
+        result = math.exp(x)
+    else:
+        # 下面改成n阶导数 ##############################################################################
+        result = math.exp(x)
+    return result
+
+
+if __name__ == '__main__':
+    ##############################################################################################################
+    list_x = [0.0,0.2,0.4,0.6,0.8]                 # 已给出的x数值，与y数值对应
+    list_y = [1.0,1.2214,1.4918,1.8221,2.2255]     # 已给出的y数值，与x数值对应
+    x_to_predict = 0.12                            # 要预测的x值
+    # 记得改下面的点数(第二个参数)，n是点数，阶数为n-1 ################################
+    print("三点牛顿前插结果为%f, 截断误差%f" % NewtonForwardInterpolation(x_to_predict, 3, list_x, list_y,FxDiff_n1))
+    print("四点牛顿前插结果为%f, 截断误差%f" % NewtonForwardInterpolation(x_to_predict, 4, list_x, list_y,FxDiff_n1))
+    print("三点牛顿基本插值结果为%f" % NewtonBaseInterpolation(x_to_predict, 3, list_x, list_y))
+    print("四点牛顿基本插值结果为%f" % NewtonBaseInterpolation(x_to_predict, 4, list_x, list_y))

按方法整理/插值-线性-抛物线-拉格朗日多项式.py

@@ -0,0 +1,79 @@
+import math
+
+# 获取与待求x最接近的两个点
+def GetClosestTwo(x,list_x):
+    for i in range(0, len(list_x)):
+        if x < list_x[i]:
+            return i-1, i
+    return len(list_x)-2, len(list_x)-1
+
+# 获取与待求x最接近的三个点
+def GetClosestThree(x,list_x):
+    if x < list_x[1]:
+        return 0, 1, 2
+    for i in range(3, len(list_x)):
+        if x < list_x[i]:
+            return i-2, i-1, i
+    return len(list_x)-3, len(list_x)-2, len(list_x)-1
+
+# 线性插值余项计算
+def LinearRegression(x,list_x,FxDiff_n):
+    i, j = GetClosestTwo(x,list_x)
+    ks = max([abs(FxDiff_n(list_x[i],2)), abs(FxDiff_n(list_x[j],2))])
+    omg = (x-list_x[i])*(x-list_x[j])
+    return abs(ks*omg/2)
+
+# 抛物线插值余项计算
+def ParabolaRegression(x,list_x,FxDiff_n):
+    i, j, k = GetClosestThree(x,list_x)
+    ks = max([abs(FxDiff_n(list_x[i],3)), abs(FxDiff_n(list_x[j],3)), abs(FxDiff_n(list_x[k],3))])
+    omg = (x-list_x[i])*(x-list_x[j])*(x-list_x[k])
+    return abs(ks*omg/6)
+
+# 线性插值
+def LinearInterpolation(x,list_x,list_y,FxDiff_n):
+    i, j = GetClosestTwo(x,list_x)
+    result = list_y[i] + (x - list_x[i]) * (list_y[j] - list_y[i]) / (list_x[j] - list_x[i])
+    r = LinearRegression(x,list_x,FxDiff_n)
+    return (result,r)
+
+# 抛物线插值
+def ParabolaInterpolation(x,list_x,list_y,FxDiff_n):
+    i, j, k = GetClosestThree(x,list_x)
+    result  = list_y[i] * (x - list_x[j]) * (x - list_x[k]) / (list_x[i] - list_x[j]) / (list_x[i] - list_x[k])
+    result += list_y[j] * (x - list_x[i]) * (x - list_x[k]) / (list_x[j] - list_x[i]) / (list_x[j] - list_x[k])
+    result += list_y[k] * (x - list_x[i]) * (x - list_x[j]) / (list_x[k] - list_x[i]) / (list_x[k] - list_x[j])
+    r = ParabolaRegression(x,list_x,FxDiff_n)
+    return (result,r)
+
+# 拉格朗日插值
+def LagrangeInterpolation(x,list_x,list_y):
+    result = 0
+    for i in range(0, len(list_x)):
+        temp = 1
+        for j in range(0, len(list_x)):
+            if i != j:
+                temp *= (x - list_x[j]) / (list_x[i] - list_x[j])
+        result += temp * list_y[i]
+    return result
+
+# 定义原函数和其导函数计算结果，用于计算插值的截断误差，如果不需要则不用管
+def FxDiff_n1(x,n):
+    result = 0
+    if n == 0:
+        # 下面改成原函数 ############################################################
+        result = math.log(x)
+    else:
+        # 下面改成n阶导数 ##############################################################################
+        result = (-1)**(n+1) * math.factorial(n-1) / (x**n)
+    return result
+
+if __name__ == "__main__":
+    ##############################################################################
+    list_x = [10,11,12,13]                    # 已给出的x数值，与y数值对应
+    list_y = [2.3026,2.3979,2.4849,2.5649]    # 已给出的y数值，与x数值对应
+    x_to_predict = 11.75                      # 要预测的x值
+    print("线性插值结果为%f, 截断误差%f" % LinearInterpolation(x_to_predict,list_x,list_y,FxDiff_n1))
+    print("抛物线插值结果为%f, 截断误差%f" % ParabolaInterpolation(x_to_predict,list_x,list_y,FxDiff_n1))
+
+    # print(LagrangeInterpolation(x_to_predict,list_x,list_y))

按方法整理/数值积分-复合梯形-复合辛普生法.py

@@ -0,0 +1,38 @@
+
+
+# 复合Newton-Cotes公式 复合牛顿-柯特斯公式
+# n等分参数，x1到x2的区间，type=1表示梯形法，type=2表示辛普生法
+def CompositeNewtonCotes(x_start,x_end,fx,n, type):
+    if type == 1:
+        h = (x_end - x_start) / n
+        result = 0
+        for i in range(n):
+            result += (fx(x_start + i * h) + fx(x_start + (i + 1) * h))
+        result *= (h / 2)
+        return result
+    elif type == 2:
+        h = (x_end - x_start) / n
+        result = -fx(x_start) + fx(x_end)
+        for i in range(n):
+            result += (4 * fx(x_start + (i + 0.5) * h) + 2 * fx(x_start + i * h))
+        result *= (h / 6)
+        return result
+    
+
+# 积分原函数 ##############################################################
+def fx(x):
+    if x == 0:
+        x = 1e-10  # Avoid division by zero  #如果x能为0，注释掉这行##############
+        pass
+    return x/(4+x**2)   #把函数改成题干的形式###################
+
+if __name__ == "__main__":
+    ##############################################################################################################
+    x_start = 3.0  # 积分下限
+    x_end = 6.0    # 积分上限
+    # 复合梯形公式，点数为n+1
+    print("复合梯形公式\n", CompositeNewtonCotes(x_start,x_end,fx,8, 1)) #8等分，1代表是梯形公式####################
+    # 复合辛普生公式，点数为2n+1
+    print("复合辛普生公式\n", CompositeNewtonCotes(x_start,x_end,fx,4, 2)) #4等分，2代表是辛普生公式###############
+
+

按方法整理/数值积分-逐次分半梯形递推公式.py

@@ -0,0 +1,38 @@
+import math
+
+# 逐次分半梯形递推公式
+def SplitTrapezoidal(a,b,fx,err):
+    count = 1
+    t1 = (b-a)*(fx(a)+fx(b))/2
+    print(f"t{count}={t1}")
+    k = 1
+    while True:
+        tmp = 0
+        for i in range(1, 2**(k-1)+1):
+            tmp += fx(a + (b-a)*(2*i-1)/(2**k))
+        t2 = t1/2+(b-a)*tmp/(2**k)
+        count *= 2
+        print(f"t{count}={t2}")
+        if abs(t2-t1) < err:
+            break
+        t1 = t2
+        k += 1
+    return t2,k
+
+
+# 积分原函数 ##############################################################
+def fx(x):
+    if x == 0:
+        x = 1e-10  # Avoid division by zero  #如果x能为0，注释掉这行##############
+        pass
+    return 1/x   #把函数改成题干的形式###################
+
+
+if __name__ == "__main__":
+    ##############################################################################################################
+    x_start = 1           # 积分下限
+    x_end = 3             # 积分上限
+    err = 1e-2            # 精度要求   P106
+
+    result,k = SplitTrapezoidal(x_start, x_end, fx, err)
+    print(f"Result: {result},k={k}")

按方法整理/数值积分-龙贝格算法.py

@@ -0,0 +1,71 @@
+import math
+
+# 龙贝格方法 积分
+def Romberg(a, b, fx, err):
+    table = [[],[],[],[]]
+    t00 = (b-a)*(fx(a)+fx(b))/2
+    table[0].append(t00)
+    print(f"t0(0)={t00}")
+    t01 = t10 = t11 = t20 = t21 = t30 = t31 = 0
+    k = 1
+    while True:
+        tmp = 0
+        for i in range(1, 2**(k-1)+1):
+            tmp += fx(a + (b-a)*(2*i-1)/(2**k))
+        t01 = t00/2+(b-a)*tmp/(2**k)
+        print(f"t0({k})={t01}")
+        table[0].append(t01)
+        
+        if k>1:
+            t11 = (4*t01-t00)/3
+            print(f"t1({k-1})={t11}")
+            table[1].append(t11)
+            if k>2:
+                t21 = (16*t11-t10)/15
+                print(f"t2({k-2})={t21}")
+                table[2].append(t21)
+                if k>3:
+                    t31 = (64*t21-t20)/63
+                    print(f"t3({k-3})={t31}")
+                    table[3].append(t31)
+                    if abs(t31-t30) < err:
+                        break
+                    t30 = t31
+                else:
+                    t30 = (64*t21-t20)/63
+                    print(f"t3(0)={t30}")
+                    table[3].append(t30)
+                t20 = t21
+            else:
+                t20 = (16*t11-t10)/15
+                print(f"t2(0)={t20}")
+                table[2].append(t20)
+            t10 = t11
+        else:
+            t10 = (4*t01-t00)/3
+            print(f"t1(0)={t10}")
+            table[1].append(t10)
+        t00 = t01
+        k += 1
+    print("Romberg table:")
+    for i in range(len(table)):
+        print(f"t{i}: {table[i]}")
+    return t31, k
+
+
+# 积分原函数 ##############################################################
+def fx(x):
+    if x == 0:
+        x = 1e-10  # Avoid division by zero  #如果x能为0，注释掉这行##############
+        pass
+    return math.sin(x)/x   #把函数改成题干的形式###################
+
+
+if __name__ == "__main__":
+    ##############################################################################################################
+    x_start = 0       # 积分下限
+    x_end = 1         # 积分上限
+    err = 0.5e-6      # 精度要求   P113
+
+    result, k = Romberg(x_start, x_end, fx, err)
+    print(f"Result: {result}, k={k}")

按方法整理/矩阵-G-S得G.py

@@ -0,0 +1,91 @@
+#模 范数
+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
+
+# 计算矩阵的逆矩阵
+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(A):",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 FixInv(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):
+            t = [row[:j] + row[j+1:] for row in (A[:i] + A[i+1:])]
+            B[j][i] = ((-1) ** (i + j)) * Det(t)
+    det = Det(A)
+    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
+
+
+
+if __name__ == "__main__":
+    ##########################################################################
+    # 把矩阵换成题干的矩阵 #########################
+    A = [
+        [1,2,-2],
+        [1,1,1],
+        [2,2,1]
+    ]
+
+    DL = [[0 for _ in range(len(A))] for __ in range(len(A))]
+    U = [[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:
+                DL[i][j] = A[i][j]
+            else:
+                U[i][j] = -A[i][j]
+                
+
+    
+    print(f"LU分解的D-L矩阵为: {DL}")
+    print(f"LU分解的U矩阵的逆矩阵为: {U}")
+    ID = FixInv(DL)
+    G = Dot(ID, U)
+    print(f"LU分解的G矩阵为: {G}")
+    

按方法整理/矩阵-SOR逐次超松弛迭代法.py

@@ -0,0 +1,62 @@
+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
+
+#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__":
+    ##########################################################################
+    #把矩阵改成题干的矩阵，b改成题干结果，err精度要求修改##########################
+    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                   # 松弛因子,题干要求    P201
+    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 '未收敛'}")

按方法整理/矩阵-列主元高斯消元法.py

@@ -0,0 +1,99 @@
+import math
+
+# 列主元高斯消元法
+def SovleRowMain(A,b):
+    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]
+
+        print(f"第{i+1}次交换，交换行{i+1}和行{row_max_index+1}，A矩阵为：")
+        for row in A:
+            print(row)
+        print(f"b向量为：{b}\n")
+
+        if abs(A[i][i]) < ks:
+            print("A矩阵奇异，无法进行高斯消元")
+            return None,None,None,None
+        for j in range(i + 1, n):
+            m = A[j][i] / A[i][i]
+            A[j][i] = m
+            for k in range(i + 1, n):
+                A[j][k] -= m * A[i][k]
+            b[j] -= m * b[i]
+            print(f"系数为{-1*m}用加号")
+        print("消元后的A矩阵：")
+        for row in A:
+            print(row)
+        print(f"消元后的b向量：{b}\n")
+
+    if abs(A[n - 1][n - 1]) < ks:
+        print("A矩阵奇异，无法进行高斯消元")
+        return None,None,None,None
+    
+    # 回代求解
+    b[n - 1] /= A[n - 1][n - 1]
+    for i in range(n - 2, -1, -1):
+        for j in range(i + 1, n):
+            b[i] -= A[i][j] * b[j]
+        b[i] /= A[i][i]
+
+    # 得到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 prettyPrintMatrix(matrix):
+    for row in matrix:
+        print(row)
+
+
+if __name__ == "__main__":
+    ##########################################################################
+    #把矩阵A和b改成题干要求的#####################################
+    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)
+    print("L:")
+    prettyPrintMatrix(L)
+    print("U:")
+    prettyPrintMatrix(U)
+    print("x:")
+    print(x)
+        
+ 

按方法整理/矩阵-平方根法.py

@@ -0,0 +1,62 @@
+import math
+
+# 获取下三角矩阵的索引
+# 列 行
+def getIndexFromDownMatrix(col, row):
+    if row > col:
+        row, col = col, row
+    return (1+col)*col//2+row
+
+
+# 平方根法求解
+def SqrtSolve(A,b):
+    n = len(b)
+    L = [0]*(n*(n+1)//2)
+    for j in range(n):
+        for i in range(j,n):
+            if i == j:
+                L[getIndexFromDownMatrix(i,j)] = A[getIndexFromDownMatrix(i,j)]
+                for k in range(j):
+                    L[getIndexFromDownMatrix(i,j)] -= L[getIndexFromDownMatrix(j,k)]**2
+                L[getIndexFromDownMatrix(i,j)] = L[getIndexFromDownMatrix(i,j)]**0.5
+            else:
+                L[getIndexFromDownMatrix(i,j)] = A[getIndexFromDownMatrix(i,j)]
+                for k in range(j):
+                    L[getIndexFromDownMatrix(i,j)] -= L[getIndexFromDownMatrix(i,k)]*L[getIndexFromDownMatrix(j,k)]
+                L[getIndexFromDownMatrix(i,j)] /= L[getIndexFromDownMatrix(j,j)]
+
+    # 打印下三角矩阵
+    print("下三角矩阵 L:")
+    for i in range(n):
+        L_row = []
+        for j in range(n):
+            if j <= i:
+                L_row.append(L[getIndexFromDownMatrix(i,j)])
+            else:
+                L_row.append(0)
+        print(L_row)
+    # print(L)
+    for i in range(n):
+        for k in range(i):
+            b[i] -= L[getIndexFromDownMatrix(i,k)]*b[k]
+        b[i] /= L[getIndexFromDownMatrix(i,i)]
+    # 打印 b 向量
+    print("y 向量:")
+    print(b)
+    for i in range(n-1,-1,-1):
+        for k in range(i+1,n):
+            b[i] -= L[getIndexFromDownMatrix(k,i)]*b[k]
+        b[i] /= L[getIndexFromDownMatrix(i,i)]
+    return b
+
+#把A,b换成题干的数值###########################################
+if __name__ == "__main__":
+    ##########################################################################
+    # 储存下三角矩阵 a11, a21, a22, a31, a32, a33 ...，按照这个格式写
+    A = [4,2,2,-2,-3,14]
+
+
+    b = [10,5,4]    # b向量，常数项
+
+    # print("x:")
+    print("x: \n",SqrtSolve(A,b))

按方法整理/矩阵-模数-范数-点积.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")))     # 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))
+    print("Norm of A with v=inf:", Norm(A, float("inf")))
+    print("Norm of A with v=1:", Norm(A, 1))
+    
+    

按方法整理/矩阵-特征值-谱半径.py

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

按方法整理/矩阵-迭代改善法.py

@@ -0,0 +1,102 @@
+import math
+
+#列主元高斯消元法
+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, fake_round_num=15):
+    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,fake_round_num)
+    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,fake_round_num)[3]
+        x0 = [x0[i] + d1[i] for i in range(len(x0))]
+        print(f"第{count+1}次迭代, r{count+1} = {r1}, d{count+1} = {d1}, x{count+2} = {x0}")
+        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__":
+    ##########################################################################
+    #把矩阵改成题干的矩阵，b改成题干结果，err精度要求修改##########################
+    A = [
+        [51,82],
+        [151/3,81]
+    ]
+    b = [235,232]
+
+    err = 1e-4    # 精度要求
+    N = 1000      # 迭代次数上限
+    fake_round_num = 4    # 模拟的四舍五入精度，根据题目情况或者凑过程修改
+
+    x = IterativeMethod(A, b, err, N, fake_round_num)[0]
+    print(f"解为: {x}")
+
+# 先求范数与逆矩阵（条件数）
+

按方法整理/矩阵-追赶法.py

@@ -0,0 +1,42 @@
+import math
+
+# 追赶法求解
+def ZGsolve(A,b):
+    n = len(b)
+    beta = [0]*n
+    for i in range(n):
+        if i == 0:
+            beta[i] = A[i][2] / A[i][1]
+        else:
+            beta[i] = A[i][2] / (A[i][1] - A[i][0]*beta[i-1])
+    print("beta:")
+    print(beta[:-1])
+    for i in range(n):
+        if i == 0:
+            b[i] = b[i] / A[i][1]
+        else:
+            b[i] = (b[i] - A[i][0]*b[i-1]) / (A[i][1] - A[i][0]*beta[i-1])
+    print("y:")
+    print(b)
+    for i in range(n-2,-1,-1):
+        b[i] = b[i] - beta[i]*b[i+1]
+
+    return b
+
+
+if __name__ == "__main__":
+    ##########################################################################
+    #把A,b换成题干的数值###########################################
+    # 储存追赶法A矩阵
+    A = [
+        [0,4,-1],
+        [-1,4,-1],
+        [-1,4,-1],
+        [-1,4,-1],
+        [-1,4,0]
+    ]
+    # A第一行前面有个0，最后一行后面有个0
+
+    b = [100,200,200,200,100]
+    print("x:")
+    print(ZGsolve(A,b))

按方法整理/矩阵-逆矩阵-条件数.py

@@ -0,0 +1,87 @@
+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
+
+# 计算矩阵的行列式
+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(A):",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 FixInv(A):
+    n = len(A)
+    if n == 2:
+        return [[A[1][1] / Det(A), -A[0][1] / Det(A)],
+                [-A[1][0] / Det(A), A[0][0] / Det(A)]]
+    B = [[0 for i in range(n)] for j in range(n)]
+    for i in range(n):
+        for j in range(n):
+            t = [row[:j] + row[j+1:] for row in (A[:i] + A[i+1:])]
+            B[j][i] = ((-1) ** (i + j)) * Det(t)
+    det = Det(A)
+    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)
+    inv_A = FixInv(A)
+    print(f"inv_A: {inv_A}, Norm(A, v): {Norm(A, v)}, Norm(inv_A, v): {Norm(inv_A, v)}")
+    # 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]
+    ]
+    # 把范数的种类数换成题干的要求，inf是无穷范数 #########################
+    print(f"矩阵A的条件数为: {Cond(A, float('inf')):.5f}")     # 1 1范数；2 2范数；float('inf') 无穷范数
+    # 把矩阵换成题干的矩阵 #########################
+    A = [
+        [1,2],
+        [3,4]
+    ]
+    # 把范数的种类数换成题干的要求，inf是无穷范数 ########################
+    print(f"矩阵A的条件数为: {Cond(A, float('inf')):.5f}")     # 1 1范数；2 2范数；float('inf') 无穷范数

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

按方法整理/非线性方程-二分法.py

@@ -0,0 +1,29 @@
+import math
+
+# 二分法求解方程的根
+def SolveByDivTwo(x1,x2,fx,err):
+    count = 1
+    
+    while abs(x2-x1) >= err:
+        x = (x1+x2)/2
+        print(f"k={count},a{count}={x1},b{count}={x2},x{count}={x},fx(x)={fx(x)}")
+        if fx(x) * fx(x1) < 0:
+            x2 = x
+        else:
+            x1 = x
+        count += 1
+
+    return (x1+x2)/2
+
+
+# 原函数换成题干的 #############################
+def fx(x):
+    return x**4-3*x+1
+
+if __name__ == "__main__":
+    ##############################################################################################################
+    x1 = 0.3         # x下限范围
+    x2 = 0.4         # x上限范围
+    err = 0.5e-5     # 精度要求       P125
+    root = SolveByDivTwo(x1, x2,fx,err)
+    print(f"Root: {root}")  

按方法整理/非线性方程-抛物线法.py

@@ -0,0 +1,66 @@
+import math
+
+
+#抛物线法解方程
+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           # 精度要求    P152
+    err2 = 1e-5           # 精度要求    P152
+    N = 100               # 最大迭代次数
+    x0 = 0.3              # 初始值1
+    x1 = 0.5              # 初始值2
+    x2 = 0.4              # 初始值3
+    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不收敛")

按方法整理/非线性方程-正割法.py

@@ -0,0 +1,38 @@
+import math
+
+
+#正割法计算方程
+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            # 根的误差限   P148
+    N0 = 100              # 最大迭代次数
+    x0 = 0.3              # 初始值1
+    x1 = 0.4              # 初始值2
+    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，无法收敛")

按方法整理/非线性方程-牛顿下山法.py

@@ -0,0 +1,50 @@
+import math
+
+def NewtonDownHillSolve(fx, dfx, x0, err1,err2, N0,min_t):
+    count = 0
+    print(f"k={count}, x0={x0}\n")
+    x1 = x0 + 1 + err1
+    while abs(x1 - x0) > err1 or abs(fx(x1)) > err2:
+        t = 1
+        if abs(dfx(x1)) < 1e-10:
+            print("导数为0，无法下山")
+            return None, 0
+        print(f"当前点: x0={x0}")
+        while t >= min_t:
+            x1 = x0 - t * fx(x0) / dfx(x0)
+            print(f"下山: t={t}, x1={x1}, fx(x{count+1})={fx(x1)}, fx(x{count})={fx(x0)}")
+            if abs(fx(x1)) < abs(fx(x0)):
+                break
+            t *= 0.5
+        if t < min_t:
+            print("达到最小t，下山失败")
+            return None, -2
+        # x1 = x0 - fx(x0) / dfx(x0)
+        count += 1
+        print(f"k={count}, x{count}={x1},x1-x0={abs(x1-x0)}\n")
+        if count > N0:
+            return None, -1
+        x0 = x1
+    print(f"收敛: x1={x1}, fx(x1)={fx(x1)}")
+    return x1, 1
+
+
+if __name__ == "__main__":
+    ##############################################################################################################
+    err1 = 1e-5      # 根的误差限   见P147
+    err2 = 1e-5      # 残量精度     见P147
+    N0 = 100         # 最大迭代次数
+    min_t = 1e-10    # 最小t值
+    x0 = 0.6         # 初始值
+    fx = lambda x: x**3 - x - 1       # 原函数
+    dfx = lambda x: 3*x**2 - 1        # 导函数
+
+    result, status = NewtonDownHillSolve(fx, dfx, x0, err1, err2, N0, min_t)
+    if status == 1:
+        print(f"收敛 解为: {result}")
+    elif status == -1:
+        print("不收敛")
+    elif status == -2:
+        print("下山失败")
+    else:
+        print("导数为0，无法收敛")

按方法整理/非线性方程-牛顿法.py

@@ -0,0 +1,50 @@
+import math
+
+
+# 牛顿方法求解方程
+def NewtonSolve(fx, dfx, x0, err, N0):
+    count = 0
+    print(f"k={count}, x0={x0}")
+    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
+        print(f"k={count}, x{count}={x1},x1-x0={abs(x1-x0)}")
+        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
+
+# 定义原函数 ###############################
+def fx(x):
+    return x**2 + 10*math.cos(x)
+
+# 定义其导函数 ###############################
+def dfx(x):
+    return 2*x - 10*math.sin(x)
+
+if __name__ == "__main__":
+    ##############################################################################################################
+    err = 1e-5             # 根的误差限   见P141
+    N0 = 100               # 最大迭代次数
+    x0 = 1.6               # 可以指定初始值
+    # x0 = FindRootZone(fx, 0, 2, 0.1) # 也可以用找根函数给定范围找根的缩小范围
+
+    result,status = NewtonSolve(fx, dfx, x0, err,N0)
+    if status == 1:
+        print(f"fx收敛 解为: {result}")
+    elif status == -1:
+        print("fx不收敛")
+    else:
+        print("fx导数为0，无法收敛")

按方法整理/非线性方程-画图找零点.py

@@ -0,0 +1,40 @@
+import math
+import matplotlib.pyplot as plt
+
+
+# 绘制函数图像 并标记可能的零点
+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=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]:.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()
+    return x, y
+
+
+
+# 把原函数换成题干的形式，自己注意定义域########################
+def fx(x):
+    return math.exp(x)-math.sin(x)
+
+
+if __name__ == "__main__":
+    ##############################################################################################################
+    a = -2*math.pi             # 左边界
+    b = math.pi                # 右边界
+    step = 0.00001             # 步长
+    print(f"边界点: {a}, {b}")
+    print(f"步长: {step}")
+    x, y = DrawGraph(a, b, step)

按方法整理/非线性方程-迭代法.py

@@ -0,0 +1,46 @@
+def fd1(x):
+    return 1+1/x**2
+
+def fd2(x):
+    return 1/(x-1)**0.5
+
+# 迭代法求解方程
+def Renew(x,fd,err):
+    count=0
+    i = 0
+    try:
+        while True:
+            xk = fd(x)
+            print(f"当前迭代值: {xk}, 上一次迭代值: {x}, 误差: {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 = xk-x
+            x = xk
+    except Exception as e:
+        print(f"发生错误: {e}")
+        return None
+    return None
+
+
+# 定义迭代公式 x = fd(x)  ###############################
+def fd(x):
+    return 1+1/x**2
+
+if __name__ == "__main__":
+    ##############################################################################################################
+    x0 = 1.5         # 初始值
+    err = 1e-5        # 精度要求          P136
+
+    result = Renew(x0, fd, err)
+    if result is not None:
+        print(f"收敛 解为: {result:.5f}")
+    else:
+        print("不收敛")
+