threesamplet.cpp

基于vc++的三次样条类函数的计算 ThreeSampleT

// ThreeSampleT.cpp: implementation of the CThreeSampleT class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "StudyThreeSampleT.h"
#include "ThreeSampleT.h"
#include "math.h"
#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CThreeSampleT::CThreeSampleT()
{

}

CThreeSampleT::~CThreeSampleT()
{

}

void CThreeSampleT::swap(double &a, double &b)
{
 double temp; 
 temp=a; 
 a=b; 
 b=temp; 
}
//求出一阶导数B
bool CThreeSampleT::GuassSolve(double *A, double *B, double *X, int n)
{
  int i,j,k,row; 
  for(k=0;k<n-1;k++) 
  {  
    double max=0.0; 
    for(i=k;i<n;i++) 
	{  
      if(fabs(A[i*n+k])>max) 
	  {  
        max=fabs(A[i*n+k]); 
        row=i; 
	  } 
	} 
    if(max==0.0) return false; 
    if(row!=k) 
	{ 
      for(j=0;j<n;j++) 
        swap(A[row*n+j],A[k*n+j]); 
      swap(B[row],B[k]); 
	} 
    for(i=k+1;i<n;i++) 
	{  
        double m=A[i*n+k]/A[k*n+k]; 
        for(j=k+1;j<n;j++) 
		{  
           A[i*n+j]=A[i*n+j]-m*A[k*n+j]; 
		}  
        B[i]=B[i]-m*B[k]; 
	}  
  }  
  X[n-1]=B[n-1]/A[n*n-1]; 
  for(k=n-2;k>=0;k--) 
  { 
     double sum=0.0; 
     for(j=k+1;j<n;j++) 
       sum+=A[k*n+j]*X[j]; 
     X[k]=(B[k]-sum)/A[k*n+k]; 
  } 
return true; 
}
//依据边界条件来重新修正参数
void CThreeSampleT::Apend(double* miu, double *rmd, double *d, double *x, double *y, double *h, double pf1, double pf2, double ppf1, double ppf2, int n, int flag)
{
  switch(flag) 
  {  
    case 0: 
      miu[n-1]=1; 
      rmd[0]=1; 
      d[0]=6*((y[1]-y[0])/(x[1]-x[0])-pf1)/h[0]; 
      d[n-1]=6*(pf2-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]))/h[n-2]; 
    break; 
    case 1: 
      rmd[0]=miu[n-1]=0; 
      d[0]=2*ppf1; 
      d[n-1]=2*ppf2; 
    break; 
    case 2: 
       rmd[0]=rmd[n-1]=h[0]/(h[n-1]+h[0]); 
       miu[0]=miu[n-1]=1-rmd[n-1]; 
       d[0]=d[n-1]=6*((y[1]-y[0])/(x[1]-x[0])-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]))/(h[0]+h[n-2]); 
    break; 
  } 
}
//x,y,分别是样条的坐标，n代表点的个数
void CThreeSampleT::T(double *x, double *y,double* M,int n)
{
int flag=0; 
double   pf1=3.0; 
double   pf2=-4.0; 
double   ppf1=0; 
double   ppf2=0; 
double * h=new double [n]; //h应该是存放区间的大小
double * miu=new double [n]; //minu存放一阶的数据
double * rmd=new double [n]; 
double * d=new double[n]; 

memset(d,0,n*sizeof(double)); 
memset(rmd,0,n*sizeof(double)); 
memset(miu,0,n*sizeof(double)); 
//存放区间的长度
for(int i=1;i<n;i++) 
{ 
	h[i-1]=x[i]-x[i-1]; 
} 
//
for(i=1;i<n-1;i++) 
{ 
	rmd[i]=h[i]/(h[i-1]+h[i]); //三次样条中求的其中一个系数
	miu[i]=h[i-1]/(h[i-1]+h[i]); //三次样条中求的其中另外一个系数
	d[i]=6*((y[i+1]-y[i])/(x[i+1]-x[i])-(y[i]-y[i-1])/(x[i]-x[i-1]))/(h[i-1]+h[i]); 
   //求出ｄ
} 
//依靠边界条件来修正前面计算的系数　　rmd，miu，d
Apend(miu,rmd,d,x,y,h,pf1,pf2,ppf1,ppf2,n,0); 
//构建一个矩阵A[n*n]
//构建一个一阶导数M
double * A=new double [n*n]; 
//,* M=new double [n]
memset(A,0,n*n*sizeof(double)); 

//给矩阵赋值
for(i=1;i<n-1;i++) 
{ 
	A[i*n+i]=2; //
	A[i*n+i-1]=miu[i]; 
	A[i*n+i+1]=rmd[i]; 
} 

if(n==2) //如果数据只有2个点时的矩阵A要特别处理
{ 
	A[0]=2; 
	A[n-2]=miu[0]; 
	A[1]=rmd[0]; 
	A[(n-1)*n+n-1]=2; 
	A[(n-1)*n+n-2]=miu[n-1]; 
	A[(n-1)*n+1]=rmd[n-1]; 
} 
else //另外对两个对角线的特殊点进行处理
{ 
    A[0]=2; 
    A[1]=rmd[0]; 
    A[(n-1)*n+n-1]=2; 
    A[(n-1)*n+n-2]=miu[n-1]; 
} 
//求出一阶导数M
GuassSolve(A,d,M,n); 
//将求出一阶导数打印出来
//for(i=0;i<n;i++) 
//{ 
//	printf("M[%d]=%f\n",i+1,M[i]); 
//} 

}
//x,y，构成样条的基本点
void CThreeSampleT::Initiate(double *x, double *y)
{
    x[0]=27.7; 
	y[0]=4.1; 

	x[1]=28; 
	y[1]=4.3; 

	x[2]=29; 
	y[2]=4.1; 

	x[3]=30; 
	y[3]=3.0; 
}
//得到利用三次样条的插值--针对单点
double CThreeSampleT::GetInterpolationValue(double *M, double *OriginX, double *OriginY, int n, double x)
{
  double* h=new double[n];
  memset(h,0,sizeof(double)*n);
  //存放区间的长度
  for(int i=1;i<n;i++) 
  {  
	  h[i-1]=OriginX[i]-OriginX[i-1]; 
  } 
  //得到x目前所处的区间号
  int RegionNo=-1;
  
  //得到需要插值的数所在区间
  RegionNo=GetXRegionNo(OriginX,n,x);

  double S=-99999.0f;
  
  int IndexNo=RegionNo-1;//应用号是区间号-1
  int j=IndexNo;
  S=M[j]*pow(OriginX[j+1]-x,3)/(6*h[j])
	  +M[j+1]*pow(x-OriginX[j],3)/(6*h[j])
	  +(OriginY[j]-M[j]*pow(h[j],2)/6)*((OriginX[j+1]-x)/h[j])
	  +(OriginY[j+1]-M[j+1]*pow(h[j],2)/6)*((x-OriginX[j])/h[j]);
  return S;


}

int CThreeSampleT::GetXRegionNo(double *OriginX, int n, double x)
{
  int RegionNo=-1;
  for(int i=1;i<n;i++)
  {
	  if(x>OriginX[i-1]&&x<OriginX[i])
	  {
         RegionNo=i;
		 break;
	  }
	  else
	  {
		  if(x==OriginX[0])
		  {
			  RegionNo=i;
			  break;
		  }
		  else
		  {
              if(x==OriginX[i])
			  {
				  RegionNo=i;
				  break;
			  }
			  else
			  {
				  if(x==OriginX[n-1])
                  {
					  RegionNo=n-1;
					  break;
				  }

			  }
		  }
	  }
  }
  return RegionNo;
}
//针对单点的案例
double CThreeSampleT::CaseTest()
{
  int n=4; 
  double * x=new double[n]; 
  double * y=new double[n]; 
  double * M=new double[n];
  Initiate(x,y); 

  T(x,y,M,n); 
  double Temp;
  Temp=GetInterpolationValue(M,x,y,n,28.5f);
  delete x;
  x=NULL;
  delete y;
  y=NULL;
  delete M;
  M=NULL;
  return Temp;
}
//针对多点的案例
void CThreeSampleT::GetInterpolationValue_MultiP(double *M, double *OriginX, double *OriginY, int n, double *OutputX, double *OutputY,int MultiPN)
{
 double x=-9999.0f;
 for(int i=0;i<MultiPN;i++)
 {
	 x=OutputX[i];
	 OutputY[i]=GetInterpolationValue(M,OriginX,OriginY,n,x);
 }
}

void CThreeSampleT::CaseTest_MultiP()
{
  int n=4; 
  double * x=new double[n]; 
  double * y=new double[n]; 
  double * M=new double[n];
  double * OutputX=new double[n];
  double * OutputY=new double[n];
  
  memset(OutputY,0,sizeof(double)*n);

  OutputX[0]=27.9;
    OutputX[1]=28.5;
	  OutputX[2]=29.2;
	    OutputX[3]=29.6;

  Initiate(x,y); 

  T(x,y,M,n); 
  double Temp;
  GetInterpolationValue_MultiP(M,x,y,n,OutputX,OutputY,n);
  delete x;
  x=NULL;
  delete y;
  y=NULL;
  delete M;
  M=NULL;
 
}