cubic.txt

三次样条插值函数的具体应用程序

/**************************************************/
//               三次样条插值函数
/**************************************************/

Cubic_spline(double RT)               
{
	
	static double g[2][N]={{112.80,63.65,37.00,22.55,14.18,9.168,6.083,4.133,2.869,2.032,1.466,1.075,0.801,0.606,0.465},
                           {-40,-30,-20,-10,0,10,20,30,40,50,60,70,80,90,100}};
	double x[N];
    double y[N];
    double h[N],a[N],b[N],A[N],B[N],m[N],s; //s表示曲线,RT表示实测热敏电阻温度
	int i,k;          
    for( i=0;i<15;i++)
       {  
          x[i]=g[1][i];
          y[i]=g[2][i];
        }
    h[0] = x[1]-x[0];       //初始化h0,a0,b0,A0,B0
    a[0] = 1;
    b[0] = 3*(y[1]-y[0])/h[0];
    A[0] = -a[0]/2;
    B[0] = b[0]/2;
 
   for(i=0;i<N;i++)        //求hi
       h[i]=x[i+1]-x[i];
 
   for(i=1;i<N-1;i++)     //求ai,bi
    {
       a[i]=h[i-1]/(h[i-1]+h[i]);  
       b[i] = 3 * ((1-a[i])*(y[i]-y[i-1])/h[i-1] + a[i]*(y[i+1]-y[i])/h[i]);
    }
   for(i=1;i<N-1;i++)    //求Ai,Bi
    {
       A[i] = -a[i]/(2+(1-a[i])*A[i-1]); 
       B[i] = (b[i]-(1-a[i])*B[i-1])/(2+(1-a[i])*A[i-1]);
     }

     m[N-1]=(b[N-1]-(1-a[N-1])*B[N-2])/(2+(1-a[N-1])*A[N-2]);     //求Mn的值
    for(i=N-2;i>=0;i--)    //求m0,m1,-----mn-1的值
           m[i] = A[i]*m[i+1] + B[i];
    for(i=0;i<N;i++)      //找出x[k]所在的区间
      {
        if((RT>=x[i])&&(RT<=x[i+1]))
//s[k]即为x[k]所在区间的三次样条插值函数,以下即为求在x[k]处的函数值 
              s=(1+2*(RT-x[i])/(x[i+1]-x[i]))*pow((RT-x[i+1])/(x[i]-x[i+1]),2)*y[i]+(1+2*(RT-x[i+1])/(x[i]-x[i+1]))*pow((RT-x[i])/(x[i+1]-x[i]),2)*y[i+1] +
                (RT-x[i])*pow((RT-x[i+1])/(x[i]-x[i+1]),2)*m[i] +(RT-x[i+1])*pow((RT-x[i])/(x[i+1]-x[i]),2)*m[i+1] ;  
       }
    return(s);
}