基于vc的三次样条类

// 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; }