（十四）三次样条插值

 

  1 #encoding=utf-8
  2 from numpy import *
  3 import  numpy as  np
  4 import matplotlib.pyplot as plt
  5 plt.close()
  6 fig=plt.figure()
  7 plt.grid(True)
  8 plt.xlabel("X")
  9 plt.ylabel("Y")
 10 plt.title(u"三次样条", fontproperties='SimHei')
 11 #三次样条
 12 #设已知n+1个点 设方程式为yi = ai*x**3+bi*x**2+ci*x+di,则需要求解4n个参数，则需要4n个条件。
 13 # 由曲线过点得2n 个条件 ，由内点处一阶导相等得n-1个条件，由内点处二阶导相等得n-1个条件 ，令端点处二阶导为0 得2个条件 a1 = 0 an = 0
 14 
 15 
 16 
 17 # #n+1个点
 18 # xi = [3,4.5,7,9]#存储x的值
 19 # n=len(xi)-1
 20 # yi = [2.5,1,2.5,0.5]#存储y 的值
 21 
 22 # 输入相关数值
 23 n=input("请输入取到样点数目：")-1
 24 xi = []
 25 yi = []
 26 for t in range(0,n+1):
 27     inputX = "请输入第"+str(t+1)+"个点 X："
 28     inputY = "请输入第"+str(t+1)+"个点 Y："
 29     xi.append(input(inputX))
 30     yi.append(input(inputY))
 31 
 32 
 33 xi_2 = []# 存储 x**2 的值
 34 xi_3 = []#存储x**3的值
 35 #向 x**2  x**3 中存入数据
 36 for i in range(0,len(xi)):
 37     xi_2.append(xi[i]**2)
 38     xi_3.append(xi[i]**3)
 39 #将方程组转化为矩阵 使 mat1 * mat2 = mat3
 40 # 可以得到行列为3n 的矩阵
 41 m1 = [[0 for a in range(4*n)] for b in range(4*n)]
 42 m3 = [[0 for a in range(1)] for b in range(4*n)]
 43 #若 mat2 内的值定为 a1,b1 ,c1,d1,a2,b2,c2,d2,a3,b3,c3,d3...... 便可以确定mat1的值
 44 #向 m1 中存入数据
 45 #定义变量 p 用于记录向矩阵插入的行数，及后续插入
 46 p=0
 47 for j in range(0,n):
 48 
 49 #从0 到n-1 的 n 个点代入
 50     m1[p][4*j]=xi_3[j]
 51     m1[p][4*j+1]=xi_2[j]
 52     m1[p][4*j+2]=xi[j]
 53     m1[p][4*j+3]=1
 54     m3[p][0]=yi[j]
 55     p=p+1
 56 #从 1 到 n 的 n 个点代入
 57     m1[p][4*j]=xi_3[j+1]
 58     m1[p][4*j+1]=xi_2[j+1]
 59     m1[p][4*j+2]=xi[j+1]
 60     m1[p][4*j+3]=1
 61     m3[p][0]=yi[j+1]
 62     p=p+1
 63 # 中间节点的一阶导 相等 2ai + bi - 2a(i+1)-b(i+1) = 0
 64 #从 1 到n-1 的 n-1 个 点处代入
 65 for k in range(1,n):
 66     #一阶导相等
 67     m1[p][4*(k-1)]= 3*xi_2[k]
 68     m1[p][4*(k-1)+1]= 2*xi[k]
 69     m1[p][4*(k-1)+2]= 1
 70 
 71     m1[p][4*k]=-3*xi_2[k]
 72     m1[p][4*k+1]=-2*xi[k]
 73     m1[p][4*k+2]=-1
 74     p=p+1
 75     #二阶导相等
 76     m1[p][4*(k-1)]=6*xi[k]
 77     m1[p][4*(k-1)+1]=2
 78 
 79     m1[p][4*k]=-6*xi[k]
 80     m1[p][4*k+1]=-2
 81     p=p+1
 82 #  代入条件 端点三阶导为零a1 = 0 an = 0
 83 m1[p][0] = 6
 84 p=p+1
 85 m1[p][4*(n-1)] = 6
 86 p=p+1
 87 #将list转化为 矩阵
 88 mat1 = np.matrix(m1)
 89 mat3 = np.matrix(m3)
 90 #求mat2
 91 _mat1 = mat1.getI()
 92 mat2=_mat1*mat3
 93 # mat2=mat3/mat1
 94 #将矩阵转化为list提取数据
 95 m2=mat2.tolist()
 96 #整理求得的曲线
 97 line=[]
 98 #用于收集区间
 99 interval=[]
100 for q in range(0,n):
101     interval.append(np.linspace(xi[q],xi[q+1]))
102     a=m2[q*4]
103     b=m2[q*4+1]
104     c=m2[q*4+2]
105     d=m2[q*4+3]
106     line.append(a*interval[q]**3+b*interval[q]**2+c*interval[q]+d)
107     plt.plot(interval[q],line[q])
108 plt.show()

 

转载于:https://www.cnblogs.com/the-wang/p/8021555.html