From f590758df3974139498f3f64ae1b5eb95d14ee88 Mon Sep 17 00:00:00 2001 From: 123 <629825095@qq.com> Date: Fri, 13 Jun 2025 22:24:52 +0800 Subject: [PATCH] 111 --- 69-5.py | 37 +++++++++++++++++++++++++++++++------ 1 file changed, 31 insertions(+), 6 deletions(-) diff --git a/69-5.py b/69-5.py index 36e56bf..d5f65f2 100644 --- a/69-5.py +++ b/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)