三次样条插值

一、样条函数的定义

样条函数属于分段光滑插值，他的基本思想是，在由两相邻节点所构成的每一个小区间内用低次多项式来逼近，并且在各结点的连接处又保证是光滑的(即导数连续)。

设在区间[a,b]上给定一组结点X:

,和一组对应的函数值。若函数S(x)满足下列条件：

(1)在每一个子区间(k=1,2,...n)上，S(x)是一个不超过三次的多项式。

(2)在每一个结点上满足S(xi)=yi,i=0,1,...,n。

(3)S(x)在区间[a,b]上为二次连续可微函数。

则称S(x)为结点X上插值与Y的三次样条插值。

二、三次样条函数的构造

在工程上，构造三次样条插值函数通常有两种方法：一是以给定插值结点处得二阶导数值作为未知数来求解，而工程上称二阶导数为弯矩，因此，这种方法成为三弯矩插值。二是以给定插值结点处得一阶导数作为未知数来求解，而一阶导数右称为斜率，因此，这种方法称为三斜率插值。

三斜率插值法

根据定义，三次样条函数在插值结点处一阶导数应存在。因此设各结点处的一阶导数为：

。利用两点埃尔米特插值公式，就可以得到样条插值函数S(x)在子区间上的表达式为：

 (1)

其中：。由此可知只要确定各结点的一阶导数值mk(k=0,1,2....n)，则各子区间上的三次样条插值函数S(x)也就确定了。

由于S(x)在区间[a,b]上的二阶导数是连续的，即在各结点的左右两子区间上的S(x)虽然不同，但在连接点的二阶导数存在，即在连接点处的二阶左导数与二阶右导数相等：。分别求子区间[xk-1,xk]右端点xk上的二阶左导数，以及子区间[xk,xk+1]左端点xk上的右导数，即可得到：

，(2)与 ，(3) 整理后可得：，k=1,2,...n-1,(4)。

其中，并且ak+bk=1。由于有n+1个插值点，因此需要确定n+1个m(m0,m1,...,mn)。而现在方程组只有n-1个方程(公式(4))，因此还需要另外补充两个条件财能唯一确定n+1个m。在实际应用中，这两个条件可以根据实际的物理意义所给出的边界条件来得到。实际问题中常用的三种边界条件如下：

(1)给定区间两个端点处的一阶导数值，即

则根据方程组(4)可得到新的方程组：

方程组的系数矩阵为，可知此矩阵为三对角矩阵。因此方程组可用追赶法求解。

解出mk后，即得到各结点的一阶导数值，将mk带入各结点的二阶导数值得表达式公式(2)或(3)可求得各结点的二阶导数值，将mk带入各区间上S(x)的表达式公式(1)即可得到各子区间上的三次样条插值函数S(x)。

(2)给定区间两个端点处得二阶到数值，即

，由此可得两个补充方程，其中。与公式(4)联立可得如下方程组:

(3)插值函数为周期函数

yn=y0，且，由此可得两个补充方程为，其中

，并且an+bn=1。最后可得方程组：

在此给出第一种边界条件下三次样条插值的算法实现：

//三次样条插值//x,y:[in] 已知点对应的函数值//N:[in] 已知点的个数//dy0:[in] x[0]的一阶导数//dyn:[in] x[n-1]的一阶导数//t:[in] 待求点的值//z:[out] t中的点所对应的函数值//nLen:[in] t与z的长度void Spline3(double *x,double *y,int N,double dy0,double dyn,double *t,double *z,int nLen){double h0,h1;double ak,bk,ck;double *pBuffer;double *dy;CBaseFunc::Allocate(pBuffer,(N+1)*N);CBaseFunc::Allocate(dy,N);memset(pBuffer,0,sizeof(double)*(N+1)*N);memset(dy,0,sizeof(double)*N);h0=x[1]-x[0];h1=x[2]-x[1];ak=h1/(h0+h1);bk=h0/(h0+h1);ck=3*(h0*(y[2]-y[1])/h1+h1*(y[1]-y[0])/h0)/(h0+h1);//求方程组的系数矩阵pBuffer[0*(N+1)+0]=1;pBuffer[0*(N+1)+N]=dy0;pBuffer[1*(N+1)+1]=2;pBuffer[1*(N+1)+2]=bk;pBuffer[1*(N+1)+N]=ck-ak*dy0;for(int i=2;i<N-2;i++){h0=h1;h1=x[i+1]-x[i];ak=h1/(h0+h1);bk=h0/(h0+h1);ck=3*(h0*(y[i+1]-y[i])/h1+h1*(y[i]-y[i-1])/h0)/(h0+h1);pBuffer[i*(N+1)+i-1]=ak;pBuffer[i*(N+1)+i]=2;pBuffer[i*(N+1)+i+1]=bk;pBuffer[i*(N+1)+N]=ck;}h0=h1;h1=x[N]-x[N-1];ak=h1/(h0+h1);bk=h0/(h0+h1);ck=3*(h0*(y[N]-y[N-1])/h1+h1*(y[N-1]-y[N-2])/h0)/(h0+h1);pBuffer[(N-2)*(N+1)+N-3]=ak;pBuffer[(N-2)*(N+1)+N-2]=bk;pBuffer[(N-2)*(N+1)+N]=ck-bk*dyn;pBuffer[(N-1)*(N+1)+N-1]=1;pBuffer[(N-1)*(N+1)+N]=dyn;//追赶法解三对角矩阵,解为各插值点的一阶导数值TridiagonalMatrix(pBuffer,N,dy);int nPos=0;memset(z,0,sizeof(double)*nLen);for(int k=0;k<nLen;k++){//寻找各待插值点的位置if(t[k]>=x[N-1]){nPos=N-2;}else{nPos=0;while(nPos+1<N){if(t[k]>=x[nPos] && t[k]<=x[nPos+1]) {break;}nPos++;}}//根据公式计算待插值点的值if(nPos+1<N){h0=x[nPos+1]-x[nPos];z[k]=dy[nPos]*pow(t[k]-x[nPos+1],2.0)*(t[k]-x[nPos])/pow(h0,2.0)+dy[nPos+1]*pow(t[k]-x[nPos],2.0)*(t[k]-x[nPos+1])/pow(h0,2.0)+y[nPos]*pow(t[k]-x[nPos+1],2.0)*(2*(t[k]-x[nPos])+h0)/pow(h0,3.0)+y[nPos+1]*pow(t[k]-x[nPos],2.0)*(-2*(t[k]-x[nPos+1])+h0)/pow(h0,3.0);}}CBaseFunc::Free(pBuffer);CBaseFunc::Free(dy);}//追赶法解三对角方程//a[k][k+1]=a[k][k+1]/a[k][k]//a[k+1][k+1]-=a[k+1][k]*a[k][k+1]bool TridiagonalMatrix(double *a,int N,double *result){for(int k=0;k<N-1;k++){a[k*(N+1)+k+1]/=a[k*(N+1)+k];a[k*(N+1)+N]/=a[k*(N+1)+k];a[(k+1)*(N+1)+k+1]-=a[(k+1)*(N+1)+k]*a[k*(N+1)+k+1];a[(k+1)*(N+1)+N]-=a[(k+1)*(N+1)+k]*a[k*(N+1)+N];}//回带求解result[N-1]=a[(N-1)*(N+1)+N]/a[(N-1)*(N+1)+N-1];for(int k=N-1-1;k>=0;k--){result[k]=(a[k*(N+1)+N]-a[k*(N+1)+k+1]*result[k+1]);}return true;}//测试函数void Spline3(){const int N=12;double x[N]={0.52,8,17.95,28.65,50.65,104.6,156.6,260.7,364.4,468,507,520};double y[N]={5.28794,13.84,20.2,24.9,31.1,36.5,36.6,31,20.9,7.8,1.5,0.2};double t[8]={4,14,30,60,130,230,450,515};double z[8];double dy0=1.86548;double dyn=-0.046115;Spline3(x,y,N,dy0,dyn,t,z,8);}
	
				
	
					   
三次样条插值
	  	   
三次样条插值
	  	   
三次样条插值
	  	   
三次样条插值
	  	   
三次样条插值
	  	   
三次样条插值
	  	   
三次样条插值
	  	   
三次样条插值
	  	   
三次样条插值
	  	   
三次样条插值 SPLINE
	  	   
二维三次样条插值
	  	   
插值法之三次样条插值
	  	   
Spline(三次样条插值)
	  	   
三次样条插值原理
	  	   
三次样条插值函数的计算
	  	   
三次样条插值函数(C++版)
	  	   
三次样条插值(Spline插值)
	  	   
Matlab 三次样条插值多项式表示
	     		  
	  	   
unity基础开发-photon服务器
	  	   
2013智能卡相关行业前景
	  	   
data masking
	  	   
jQuery中的事件绑定和去除事件绑定
	  	   
css3
	  	   
三次样条插值
	  	   
PKU 3067 Japan 树状数组
	  	   
慢慢
	  	   
JQuery validate验证消息弹出对话框
	  	   
iPad、iPhone界面元素设计标准（简洁版）
	  	   
头文件中定义 const 全局变量应注意的问题
	  	   
ExtJS事件处理
	  	   
yii-程序- 使用事务
	  	   
数据库触发器