三次样条曲线

#include<iostream>
#include<iomanip>
using namespace std;

const int M = 16;
double dKnowX[M] =
{
 -10,-8,-6,-4,-2,
 0,2,4,6,8,
 10,12,14,16,18,
 20
};
double dKnowY[M] =
{
 -715,-331,-115,-19,5,
 5,29,125,341,725,
 1325,2189,3365,4901,6845,
 9245
};
double dKnowDy[M];
double dKnowDdy[M];

const int N = 15;
double dInsertX[N] =
{
 -9,-7,-5,-3,-1,
 1,3,5,7,9,
 11,13,15,17,19
};
double dInsertY[N];
double dInsertDy[N];
double dInsertDdy[N];

 

/*
out_dY   所求函数值
out_dDy  所求插值点的一阶导数
out_dDdy 所求插值点的二阶导数
in_dX    所求插值点的插值点的X轴坐标(升序排列)
in_dKnowX,in_dKnowY  已知点的节点坐标(升序排列)
in_dKnowDy,  in_dKnowDdy 已知点的节点的一阶导和二阶导(升序排列)
in_iN 已知节点个数
*/
void CubicSpline(double out_dY[], double out_dDy[], double out_dDdy[],
     double in_dX[], int in_iM,
     double in_dKnowX[], double in_dKnowY[],
     double in_dKnowDy[], double in_dKnowDdy[],
     int in_iN)  
{
 //三次样条函数

      int   i,j;  
      double   h0,h1,alpha,beta,g,s[400];  

      in_dKnowDy[0]=-0.5;  
      h0=in_dKnowX[1]-in_dKnowX[0];  
      s[0]=3.0*(in_dKnowY[1]-in_dKnowY[0])/(2.0*h0)-in_dKnowDdy[0]*h0/4.0;  
      for(j=1;j<=in_iN-2;j++)  
   {
    h1=in_dKnowX[j+1]-in_dKnowX[j];  
    alpha=h0/(h0+h1);  
    beta=(1.0-alpha)*(in_dKnowY[j]-in_dKnowY[j-1])/h0;  
    beta=3.0*(beta+alpha*(in_dKnowY[j+1]-in_dKnowY[j])/h1);
    in_dKnowDy[j]=-alpha/(2.0+(1.0-alpha)*in_dKnowDy[j-1]);  
    s[j]=(beta-(1.0-alpha)*s[j-1]);  
    s[j]=s[j]/(2.0+(1.0-alpha)*in_dKnowDy[j-1]);  
    h0=h1;  
   }

      in_dKnowDy[in_iN-1]=(3.0*(in_dKnowY[in_iN-1]-in_dKnowY[in_iN-2])/h1+in_dKnowDdy[in_iN-1]*h1/2.0-s[in_iN-2])/(2.0+in_dKnowDy[in_iN-2]);  
      for(j=in_iN-2;j>=0;j--)
    in_dKnowDy[j]=in_dKnowDy[j]*in_dKnowDy[j+1]+s[j];  
      for(j=0;j<=in_iN-2;j++)  
    s[j]=in_dKnowX[j+1]-in_dKnowX[j];  
      for(j=0;j<=in_iN-2;j++)
   {
    h1=s[j]*s[j];
    in_dKnowDdy[j]=6.0*(in_dKnowY[j+1]-in_dKnowY[j])/h1-2.0*(2.0*in_dKnowDy[j]+in_dKnowDy[j+1])/s[j];
   }  
      h1=s[in_iN-2]*s[in_iN-2];  
      in_dKnowDdy[in_iN-1]=6.*(in_dKnowY[in_iN-2]-in_dKnowY[in_iN-1])/h1+2.*(2.*in_dKnowDy[in_iN-1]+in_dKnowDy[in_iN-2])/s[in_iN-2];  
      g=0.0;  
      for(i=0;i<=in_iN-2;i++)
   {
    h1=0.5*s[i]*(in_dKnowY[i]+in_dKnowY[i+1]);
    h1=h1-s[i]*s[i]*s[i]*(in_dKnowDdy[i]+in_dKnowDdy[i+1])/24.0;
    g=g+h1;
   }  

      for(j=0;j<=in_iM-1;j++)
   {
    if(in_dX[j]>=in_dKnowX[in_iN-1])  
     i=in_iN-2;
    else
    {
     i=0;
     while(in_dX[j]>in_dKnowX[i+1])  
      i=i+1;
    }  
          h1=(in_dKnowX[i+1]-in_dX[j])/s[i];  
          h0=h1*h1;  
          out_dY[j]=(3.0*h0-2.0*h0*h1)*in_dKnowY[i];  
          out_dY[j]=out_dY[j]+s[i]*(h0-h0*h1)*in_dKnowDy[i];  
          out_dDy[j]=6.0*(h0-h1)*in_dKnowY[i]/s[i];  
          out_dDy[j]=out_dDy[j]+(3.0*h0-2.0*h1)*in_dKnowDy[i];  
          out_dDdy[j]=(6.0-12.0*h1)*in_dKnowY[i]/(s[i]*s[i]);  
          out_dDdy[j]=out_dDdy[j]+(2.0-6.0*h1)*in_dKnowDy[i]/s[i];  
          h1=(in_dX[j]-in_dKnowX[i])/s[i];  
          h0=h1*h1;  
          out_dY[j]=out_dY[j]+(3.0*h0-2.0*h0*h1)*in_dKnowY[i+1];  
          out_dY[j]=out_dY[j]-s[i]*(h0-h0*h1)*in_dKnowDy[i+1];  
          out_dDy[j]=out_dDy[j]-6.0*(h0-h1)*in_dKnowY[i+1]/s[i];  
          out_dDy[j]=out_dDy[j]+(3.0*h0-2.0*h1)*in_dKnowDy[i+1];  
          out_dDdy[j]=out_dDdy[j]+(6.0-12.0*h1)*in_dKnowY[i+1]/(s[i]*s[i]);  
          out_dDdy[j]=out_dDdy[j]-(2.0-6.0*h1)*in_dKnowDy[i+1]/s[i];
   }  
} 

void CalculateDyDdy(double out_dDy[], double out_dDdy[], 
     double in_dPx[], double in_dPy[], int in_iN,
     double in_dBdy = 0, double in_dEdy = 0,
     double in_dBDdy = 0, double in_dEDdy = 0)
{//以估计的方法计算导数
 out_dDy[0] = in_dBdy;
 out_dDdy[0] = in_dBDdy;
 out_dDy[in_iN-1] = in_dEdy;
 out_dDdy[in_iN-1] = in_dEDdy;
 for(int i=1; i<in_iN-1; i++)
 {
  double dx = in_dPx[i+1] - in_dPx[i];
  if( 0 != dx)
  {
   out_dDy[i] = (in_dPy[i+1] - in_dPy[i]) / (2*dx);
   out_dDdy[i] = (in_dPy[i+1] + in_dPy[i-1] - 2*in_dPy[i]) / dx / dx;
  }
 }
}

void main()
{
 CalculateDyDdy(dKnowDy, dKnowDdy, dKnowX, dKnowY, M);
 CubicSpline(dInsertY, dInsertDy, dInsertDdy, dInsertX, N,
  dKnowX, dKnowY, dKnowDy, dKnowDdy, M);
 for(int i=0; i<N; i++)
 {
  cout<<dInsertX[i]<<",   "<<dInsertY[i]<<endl;
 }
 system("pause");
}