spline（三次样条插值（自然边界条件））.c

曲线拟合程序 多项式相关系数的计算方法(多项式形式1) 多项式相关系数的计算方法(多项式形式2) 最小二乘法曲线拟合 三次样条插值（自然边界条件）

//=====================================================================================
//函数说明
//函数名称：Spline
//函数功能：三次样条插值（自然边界条件）算法
//使用方法：double *x ---- 结点横坐标
//          double *y ---- 结点纵坐标
//			double *dy --- 过程变量（调用该函数之前只需为其开辟空间），个数是结点个数
//			int n ------- 结点个数
//
//          double *t ---- 插入点横坐标
//			double *z ---- 插入点拟合值，为输出值
//          double *s ---- 过程变量（调用该函数之前只需为其开辟空间），个数是插入点个数
//			int m ------- 插入点个数
//=====================================================================================
void Spline(double *x, double *y, double *dy, int n, double *t, double *z, double *s, int m)
{
	int  i, j;
	double  h0, h1, alpha, beta;
	
	//-------------------------------------------------------------------
	//计算aj和bj
	//-------------------------------------------------------------------
	dy[0] = -0.5;								  //a0
	h0 = x[1]-x[0];
	s[0] = 3.0 *(y[1]-y[0])/(2.0 *h0);  		  //b0
	
	for (j = 1; j <= n-2; j++)
	{
		h1 = x[j+1]-x[j];
		alpha = h0/(h0+h1);
		beta = (1.0-alpha)*(y[j]-y[j-1])/h0;
		beta = 3.0 *(beta+alpha *(y[j+1]-y[j])/h1);
	
		dy[j] = -alpha/(2.0+(1.0-alpha) *dy[j-1]);//aj
	
		s[j] = (beta-(1.0-alpha) *s[j-1]);
		s[j] = s[j]/(2.0+(1.0-alpha) *dy[j-1]);   //bj
		h0 = h1;
	}
	
	//-------------------------------------------------------------------
	//利用各结点上的函数值和一阶导数值计算插值点t的函数
	//-------------------------------------------------------------------
	dy[n-1] = (3.0 *(y[n-1]-y[n-2])/h1-s[n-2])/(2.0+dy[n-2]);
	
	for (j = n-2; j >= 0; j--)  	//计算结点上的一阶导数值
	{
		dy[j] = dy[j] *dy[j+1]+s[j];
	}
	
	for (j = 0; j <= n-2; j++)  	//计算hj
	{
		s[j] = x[j+1]-x[j]; 
	}
		
	for (j = 0; j <= m-1; j++)   	//计算插值点t的函数值
	{
		if (t[j] >= x[n-1])			//判断t[j]的所在区间
		{
			i = n-2;
		}
		else
		{
			i = 0;
			while (t[j] > x[i+1])
			{
				i = i+1;
			}
		}
	
		h1 = (x[i+1]-t[j])/s[i];	//求解
		h0 = h1 * h1;
		z[j] = (3.0 *h0-2.0 * h0 * h1) *y[i];
		z[j] = z[j]+s[i]*(h0-h0 * h1) *dy[i];	
		h1 = (t[j]-x[i])/s[i];
		h0 = h1 * h1;
		z[j] = z[j]+(3.0 *h0-2.0 * h0 * h1) *y[i+1];
		z[j] = z[j]-s[i]*(h0-h0 * h1) *dy[i+1];  
	}
}