111

69-5.py

@@ -1,5 +1,3 @@
-x = [0,1,2,3]
-y = [0,0,0,0]
 
 # 追赶法
 def ZGsolve(A,b):
@@ -39,15 +37,33 @@ def GetDQList(list_x, list_y):
 # 三次样条插值
 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))
@@ -56,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:  # 二阶导数边界条件
@@ -66,9 +85,13 @@ 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
 
 # 打印矩阵
@@ -80,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)