devirate.h

利用样条插值求取一组数据点在各点的导数 (只提供函数)

/////////////利用三次样条插值求取离散点的导数
struct point
{ double x,y,k,M;};
void Devirate(struct point pts[],int dt)
{
	int i,j;
	int icounter=0;
  	double a=0.00;
	double b;
	double h0;
	double h1;
	double tempa[NumberOfPt],tempb[NumberOfPt];
	pts[0].k= Five(pts[0].y,pts[1].y,pts[2].y,pts[3].y,pts[4].y,
				   pts[0].x,pts[1].x,pts[2].x,pts[3].x,pts[4].x);//计算首节点斜率
 
	pts[dt-1].k = Five(pts[dt-1].y,pts[dt-2].y,pts[dt-3].y,pts[dt-4].y,pts[dt-5].y,
					  pts[dt-1].x,pts[dt-2].x,pts[dt-3].x,pts[dt-4].x,pts[dt-5].x);//计算尾节点
	a=0.00;
	b=pts[0].y;
	h0=pts[1].x-pts[0].x;
   	tempa[0] = 0.00;
	tempb[0] = pts[0].k;
		
	 ////////////////应用三转角方程求型值点的一阶导////////////////////////
	for(j = 1; j<dt-1; j++)
	{
		h1=pts[j+1].x-pts[j].x;
		a=h0/(h0+h1) ;
		b=3.0*((1.0-a)*(pts[j].y-pts[j-1].y)/h0+a*(pts[j+1].y-pts[j].y)/h1) ;
		tempa[j] = -a/(2.0+(1.0-a)*tempa[j-1]) ;
		tempb[j]=(b-(1.0-a)*tempb[j-1])/(2.0+(1.0-a)*tempa[j-1]);
		h0=h1;
	}

	for(i = dt-2;i>0;i--)
	{
		pts[i].k=tempa[i]*pts[i+1].k+tempb[i];
	}
	 //////////////////应用三转角方程求型值点的二阶导/////////////////////////////////
	  ///应用三转角方程求型值点的二阶导
	pts[0].M= Five(pts[0].k,pts[1].k,pts[2].k,pts[3].k,pts[4].k,
				pts[0].x,pts[1].x,pts[2].x,pts[3].x,pts[4].x);//计算首节点斜率
	pts[dt-1].M = Five(pts[dt-1].k,pts[dt-2].k,pts[dt-3].k,pts[dt-4].k,pts[dt-5].k,
				pts[dt-1].x,pts[dt-2].x,pts[dt-3].x,pts[dt-4].x,pts[dt-5].x); //计算尾节点
	a = 0.00;
	b = pts[0].k;
	h0 = pts[1].x - pts[0].x;
	tempa[0] = 0.00;
	tempb[0] = pts[0].M;

	for(j = 1; j<dt-1; j++)
	{
		h1=pts[j+1].x-pts[j].x;
		a=h0/(h0+h1) ;
		b=3.0*((1.0-a)*(pts[j].k-pts[j-1].k)/h0+a*(pts[j+1].k-pts[j].k)/h1) ;
		tempa[j] = -a/(2.0+(1.0-a)*tempa[j-1]) ;
		tempb[j]=(b-(1.0-a)*tempb[j-1])/(2.0+(1.0-a)*tempa[j-1]);
		h0=h1;
	}
	for(i = dt-2;i>0;i--)
	{
		pts[i].M=tempa[i]*pts[i+1].M+tempb[i];
	} 
 }


//////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////五点求导///////////////////////////////////////////////

double Five(double f0, double f1,double f2, double f3,double f4,double x0, double x1,double x,double x3,double x4)
{
	double temp = 0.000;
	temp = ((2.0*x0-x1-x)*(x0-x3)*(x0-x4)
			+(2.0*x0-x3-x4)*(x0-x1)*(x0-x))*f0/((x0-x1)*(x0-x)*(x0-x3)*(x0-x4)) 
			+((x0-x1)*(x0-x3)*(x0-x4))*f2/((x-x0)*(x-x1)*(x-x3)*(x-x4))
			+((x0-x)*(x0-x3)*(x0-x4))*f1/((x1-x0)*(x1-x)*(x1-x3)*(x1-x4)) 
			+(x0-x1)*(x0-x)*(x0-x4)*f3/((x3-x0)*(x3-x)*(x3-x1)*(x3-x4)) 
			+(x0-x1)*(x0-x)*(x0-x3)*f4/((x4-x0)*(x4-x)*(x4-x3)*(x4-x1));
	return temp;
}