百度百科定义
插值：在离散数据的基础上插补连续函数，使得这条连续曲线经过全部离散点，同时也可以估计出函数在其他点的近似值。
样条插值：一种以 可变样条 来作出一条经过一系列点的光滑曲线的数学方法。插值样条是由一些多项式组成的，每一个多项式都是由相邻的两个数据点决定的，这样，任意的两个相邻的多项式以及它们的导数在连接点处都是连续的。
样条插值法
简单理解，就是每两个点之间确定一个函数，这个函数就是一个样条，函数不同，样条就不同，所以定义中说 可变样条，然后把所有样条分段结合成一个函数，就是最终的插值函数。
思路1 - 线性样条
两点确定一条直线，我们可以在每两点间画一条直线，就可以把所有点连起来。
显然曲线不够光滑，究其原因是因为连接点处导数不相同。
思路2 - 二次样条
直线不行，用曲线代替，二次函数是最简单的曲线。
假设4个点，x0，x1，x2，x3，有3个区间，需要3个二次样条，每个二次样条为 ax^2+bx+c，故总计9个未知数。
1. x0，x3两个端点都有一个二次函数经过，可确定2个方程
2. x1，x2两个中间点都有两个二次函数经过，可确定4个方程
3. 中间点处必须连续，需要保证左右二次函数一阶导相等
2*a1*x1+b1=2*a2*x1+b2
2*a2*x2+b2=2*a3*x2+b3
可确定2个方程，此时有了8个方程。
4. 这里假设第一方程的二阶导为0，即 a1=0，又是一个方程，共计9个方程。 【见补充】
联立即可求解。
python 实现
# encoding:utf-8import numpy as npimport matplotlib.pyplot as pltfrom pylab import mpl"""二次样条实现"""x = [3, 4.5, 7, 9]y = [2.5, 1, 2.5, 0.5]def calculateEquationParameters(x): #parameter为二维数组，用来存放参数，sizeOfInterval是用来存放区间的个数 parameter = [] sizeOfInterval=len(x)-1 i = 1 #首先输入方程两边相邻节点处函数值相等的方程为2n-2个方程 while i < len(x)-1: data = init(sizeOfInterval*3) data[(i-1)*3]=x[i]*x[i] data[(i-1)*3+1]=x[i] data[(i-1)*3+2]=1 data1 =init(sizeOfInterval*3) data1[i * 3] = x[i] * x[i] data1[i * 3 + 1] = x[i] data1[i * 3 + 2] = 1 temp=data[1:] parameter.append(temp) temp=data1[1:] parameter.append(temp) i += 1 #输入端点处的函数值。为两个方程,加上前面的2n-2个方程，一共2n个方程 data = init(sizeOfInterval*3-1) data[0] = x[0] data[1] = 1 parameter.append(data) data = init(sizeOfInterval *3) data[(sizeOfInterval-1)*3+0] = x[-1] * x[-1] data[(sizeOfInterval-1)*3+1] = x[-1] data[(sizeOfInterval-1)*3+2] = 1 temp=data[1:] parameter.append(temp) #端点函数值相等为n-1个方程。加上前面的方程为3n-1个方程,最后一个方程为a1=0总共为3n个方程 i=1 while i < len(x) - 1: data = init(sizeOfInterval * 3) data[(i - 1) * 3] =2*x[i] data[(i - 1) * 3 + 1] =1 data[i*3]=-2*x[i] data[i*3+1]=-1 temp=data[1:] parameter.append(temp) i += 1 return parameter"""对一个size大小的元组初始化为0"""def init(size): j = 0 data = [] while j < size: data.append(0) j += 1 return data"""功能：计算样条函数的系数。参数：parametes为方程的系数，y为要插值函数的因变量。返回值：二次插值函数的系数。"""def solutionOfEquation(parametes,y): sizeOfInterval = len(x) - 1 result = init(sizeOfInterval*3-1) i=1 while i<sizeOfInterval: result[(i-1)*2]=y[i] result[(i-1)*2+1]=y[i] i+=1 result[(sizeOfInterval-1)*2]=y[0] result[(sizeOfInterval-1)*2+1]=y[-1] a = np.array(calculateEquationParameters(x)) b = np.array(result) return np.linalg.solve(a,b)"""功能：根据所给参数，计算二次函数的函数值：参数:parameters为二次函数的系数，x为自变量返回值：为函数的因变量"""def calculate(paremeters,x): result=[] for data_x in x: result.append(paremeters[0]*data_x*data_x+paremeters[1]*data_x+paremeters[2]) return result"""功能：将函数绘制成图像参数：data_x,data_y为离散的点.new_data_x,new_data_y为由拉格朗日插值函数计算的值。x为函数的预测值。返回值：空"""def Draw(data_x,data_y,new_data_x,new_data_y): plt.plot(new_data_x, new_data_y, label=u"拟合曲线", color="black") plt.scatter(data_x,data_y, label=u"离散数据",color="red") mpl.rcParams['font.sans-serif'] = ['SimHei'] mpl.rcParams['axes.unicode_minus'] = False plt.title(u"二次样条函数") plt.legend(loc="upper left") plt.show()result=solutionOfEquation(calculateEquationParameters(x),y)new_data_x1=np.arange(3, 4.5, 0.1)new_data_y1=calculate([0,result[0],result[1]],new_data_x1)new_data_x2=np.arange(4.5, 7, 0.1)new_data_y2=calculate([result[2],result[3],result[4]],new_data_x2)new_data_x3=np.arange(7, 9.5, 0.1)new_data_y3=calculate([result[5],result[6],result[7]],new_data_x3)new_data_x=[]new_data_y=[]new_data_x.extend(new_data_x1)new_data_x.extend(new_data_x2)new_data_x.extend(new_data_x3)new_data_y.extend(new_data_y1)new_data_y.extend(new_data_y2)new_data_y.extend(new_data_y3)Draw(x,y,new_data_x,new_data_y)
可以看到 y 是多段二次函数拼接而成。
输出
二次样条插值连续光滑，看起来效果还行。
只是前两个点之间是条直线，因为假设a1=0，二次函数变成b1x+c1，显然是直线；
而且最后两个点之间过于陡峭 。
思路3 - 三次样条
二次函数最高项系数为0，导致变成直线，那三次函数最高项系数为0，还是曲线，插值效果应该更好。
三次样条思路与二次样条基本相同，
同样假设4个点，x0，x1，x2，x3，有3个区间，需要3个三次样条，每个三次样条为 ax^3+bx^2+cx+d，故总计12个未知数。
1. 内部节点处的函数值应该相等，这里一共是4个方程。
2. 函数的第一个端点和最后一个端点，应该分别在第一个方程和最后一个方程中。这里是2个方程。
3. 两个函数在节点处的一阶导数应该相等。这里是两个方程。
4. 两个函数在节点处的二阶导数应该相等，这里是两个方程。 【见补充】
5. 假设端点处的二阶导数为零，这里是两个方程。 【见补充】
a1=0
b1=0
python 实现
# encoding:utf-8import numpy as npimport matplotlib.pyplot as pltfrom pylab import mpl"""三次样条实现"""x = [3, 4.5, 7, 9]y = [2.5, 1, 2.5, 0.5]def calculateEquationParameters(x): #parameter为二维数组，用来存放参数，sizeOfInterval是用来存放区间的个数 parameter = [] sizeOfInterval=len(x)-1; i = 1 #首先输入方程两边相邻节点处函数值相等的方程为2n-2个方程 while i < len(x)-1: data = init(sizeOfInterval*4) data[(i-1)*4] = x[i]*x[i]*x[i] data[(i-1)*4+1] = x[i]*x[i] data[(i-1)*4+2] = x[i] data[(i-1)*4+3] = 1 data1 =init(sizeOfInterval*4) data1[i*4] =x[i]*x[i]*x[i] data1[i*4+1] =x[i]*x[i] data1[i*4+2] =x[i] data1[i*4+3] = 1 temp = data[2:] parameter.append(temp) temp = data1[2:] parameter.append(temp) i += 1 # 输入端点处的函数值。为两个方程, 加上前面的2n - 2个方程，一共2n个方程 data = init(sizeOfInterval * 4 - 2) data[0] = x[0] data[1] = 1 parameter.append(data) data = init(sizeOfInterval * 4) data[(sizeOfInterval - 1) * 4 ] = x[-1] * x[-1] * x[-1] data[(sizeOfInterval - 1) * 4 + 1] = x[-1] * x[-1] data[(sizeOfInterval - 1) * 4 + 2] = x[-1] data[(sizeOfInterval - 1) * 4 + 3] = 1 temp = data[2:] parameter.append(temp) # 端点函数一阶导数值相等为n-1个方程。加上前面的方程为3n-1个方程。 i=1 while i < sizeOfInterval: data = init(sizeOfInterval * 4) data[(i - 1) * 4] = 3 * x[i] * x[i] data[(i - 1) * 4 + 1] = 2 * x[i] data[(i - 1) * 4 + 2] = 1 data[i * 4] = -3 * x[i] * x[i] data[i * 4 + 1] = -2 * x[i] data[i * 4 + 2] = -1 temp = data[2:] parameter.append(temp) i += 1 # 端点函数二阶导数值相等为n-1个方程。加上前面的方程为4n-2个方程。且端点处的函数值的二阶导数为零，为两个方程。总共为4n个方程。 i = 1 while i < len(x) - 1: data = init(sizeOfInterval * 4) data[(i - 1) * 4] = 6 * x[i] data[(i - 1) * 4 + 1] = 2 data[i * 4] = -6 * x[i] data[i * 4 + 1] = -2 temp = data[2:] parameter.append(temp) i += 1 return parameter"""对一个size大小的元组初始化为0"""def init(size): j = 0 data = [] while j < size: data.append(0) j += 1 return data"""功能：计算样条函数的系数。参数：parametes为方程的系数，y为要插值函数的因变量。返回值：三次插值函数的系数。"""def solutionOfEquation(parametes,y): sizeOfInterval = len(x) - 1 result = init(sizeOfInterval*4-2) i=1 while i<sizeOfInterval: result[(i-1)*2]=y[i] result[(i-1)*2+1]=y[i] i+=1 result[(sizeOfInterval-1)*2]=y[0] result[(sizeOfInterval-1)*2+1]=y[-1] a = np.array(calculateEquationParameters(x)) b = np.array(result) for data_x in b: print(data_x) return np.linalg.solve(a,b)"""功能：根据所给参数，计算三次函数的函数值：参数:parameters为二次函数的系数，x为自变量返回值：为函数的因变量"""def calculate(paremeters,x): result=[] for data_x in x: result.append(paremeters[0]*data_x*data_x*data_x+paremeters[1]*data_x*data_x+paremeters[2]*data_x+paremeters[3]) return result"""功能：将函数绘制成图像参数：data_x,data_y为离散的点.new_data_x,new_data_y为由拉格朗日插值函数计算的值。x为函数的预测值。返回值：空"""def Draw(data_x,data_y,new_data_x,new_data_y): plt.plot(new_data_x, new_data_y, label=u"拟合曲线", color="black") plt.scatter(data_x,data_y, label=u"离散数据",color="red") mpl.rcParams['font.sans-serif'] = ['SimHei'] mpl.rcParams['axes.unicode_minus'] = False plt.title(u"三次样条函数") plt.legend(loc="upper left") plt.show()result=solutionOfEquation(calculateEquationParameters(x),y)new_data_x1=np.arange(3, 4.5, 0.1)new_data_y1=calculate([0,0,result[0],result[1]],new_data_x1)new_data_x2=np.arange(4.5, 7, 0.1)new_data_y2=calculate([result[2],result[3],result[4],result[5]],new_data_x2)new_data_x3=np.arange(7, 9.5, 0.1)new_data_y3=calculate([result[6],result[7],result[8],result[9]],new_data_x3)new_data_x=[]new_data_y=[]new_data_x.extend(new_data_x1)new_data_x.extend(new_data_x2)new_data_x.extend(new_data_x3)new_data_y.extend(new_data_y1)new_data_y.extend(new_data_y2)new_data_y.extend(new_data_y3)Draw(x,y,new_data_x,new_data_y)
输出
补充
微分连续性
s 代表三次样条，s'是一阶导，s''是二阶导
端点条件
上面我们对端点处的样条进行了假设，为什么呢？其实端点可以有多种不同的限制，常见有3种。
自由边界 Natural
首尾两端没有受到任何使他们弯曲的力，二次样条就是 s’=0，三次样条就是 s''=0
固定边界 Clamped
首尾两端点的微分值被指定
非节点边界 Not-A-Knot
把端点当做中间点处理，三次函数不做假设，即
不同边界的比较
参考资料：
https://blog.csdn.net/flyingleo1981/article/details/53008931 三次样条插值(Cubic Spline Interpolation)及代码实现(C语言)
https://blog.csdn.net/deramer1/article/details/79034201 三次样条插值法
---------------------
作者：努力的孔子
来源：CNBLOGS
原文：https://www.cnblogs.com/yanshw/p/11194058.html
版权声明：本文为作者原创文章，转载请附上博文链接！
内容解析By：CSDN,CNBLOG博客文章一键转载插件 标签：方程,函数,样条,插值,博客园,temp,sizeOfInterval,data From： https://www.cnblogs.com/excellentliu/p/18374680