trispline.cpp

三次样条插值程序。编译运行本程序后

//三次样条插值
#include <iostream.h>
#include <math.h>

double fun_func(double x)                               //x结点精确值计算
{
	double y;
	y=exp(-2*x);
	return y;
}

int judge(double *x,double v,int n)                    //判断所输入x值所属区间，返回i值
{
     for(int i=0;i<=n;i++)
	{
		if(x[i]<=v&&x[i+1]>v)
		{
			return i;
			cout<<i<<endl; 
		}	
	}
	 return 0;//
}
double Result_TriSpline(double *y, double *x,double *h,double *dy,double v,int i,int n)
                                                      //修改历史：int v--〉double v
			                                          //计算Y(v)的三次样条插值
{
	  double * alpha = new double[n];                 //开辟指针空间
	  double * beta = new double[n];                  //开辟指针空间
	  double * a = new double[n];                     //开辟指针空间
	  double * b = new double[n];                     //开辟指针空间
	   dy[0]=-2*exp(-2*x[0]);
	   //cout<<dy[0]<<"test"<<endl;                   //For test
       dy[n]=-2*exp(-2*x[n]);
	   //cout<<dy[n]<<"test"<<endl;                   //For test
	   a[0]=0;
	   b[0]=dy[0];
       h[0]=x[1]-x[0];
	   alpha[0]=0;                                    //修改历史：添加
	   x[n]=6;                                  //修改历史：7——〉6+6/n
	   y[n+1]=exp(-2*n);                          //修改历史：添加

	   for(int j=1;j<=n-1;j++)
	   {
		   h[j]=x[j+1]-x[j];
		   //cout<<h[j]<<"   "<<"test"<<endl;         //for test
		   alpha[j]=(h[j-1])/(h[j-1]+h[j]);
           //cout<<alpha[j]<<"   "<<"test"<<endl;     //for test
           beta[j]=3*(1-alpha[j])*(y[j]-y[j-1])/h[j-1]+3*alpha[j]*(y[j+1]-y[j])/h[j];
		   //cout<<beta[j]<<"   "<<"test"<<endl;       //for test
		   a[j]=-alpha[j]/(2+(1-alpha[j])*a[j-1]);
		   //cout<<a[j]<<"   "<<"test"<<endl;          //for test
		   b[j]=(beta[j]-(1-alpha[j])*b[j-1])/(2+(1-alpha[j])*a[j-1]);
		   //cout<<b[j]<<"   "<<"test"<<endl;          //for test		   
	   }
	
       for(int k=n-1;k>=1;k--)//
	   {
           dy[k]=a[k]*dy[k+1]+b[k];//
		   //cout<<dy[k]<<"   "<<"test"<<endl;         //  for test//
	   }
		
	double R1=(3*(x[i+1]-v)*(x[i+1]-v))/(h[i]*h[i])-(2*(x[i+1]-v)*(x[i+1]-v)*(x[i+1]-v))/(h[i]*h[i]*h[i]);
	//cout<<R1<<endl;                                   //test
	double R2=(3*(v-x[i])*(v-x[i]))/(h[i]*h[i])-(2*(v-x[i])*(v-x[i])*(v-x[i]))/(h[i]*h[i]*h[i]);//double ,^不可用~~
	//cout<<R2<<endl;//test
	double R3=((x[i+1]-v)*(x[i+1]-v))/(h[i]*h[i])-((x[i+1]-v)*(x[i+1]-v)*(x[i+1]-v))/(h[i]*h[i]*h[i]);
	//cout<<R3<<endl;//test
	double R4=((v-x[i])*(v-x[i]))/(h[i]*h[i])-((v-x[i])*(v-x[i])*(v-x[i]))/(h[i]*h[i]*h[i]);
	//cout<<R4<<endl;//test
	double R=R1*y[i]+R2*y[i+1]+R3*h[i]*dy[i]+R4*h[i]*dy[i+1];
	//cout<<i<<"  "<<y[i+1]<<"  "<<h[i]<<"  "<<dy[i]<<endl;//for  test
	
    cout<<R<<endl;                                          //for  test
	
	delete []a;
	delete []b;
	delete []alpha;
	delete []beta;
	
	return R;
}
void main()
{
	int ToF_goon=1;
	while(ToF_goon==1)
	{ 
		int n=0;double v=-1;
       while(n<=1)
	  {
	    cout<<endl<<"input the number of intervals:";
	    cin>>n;
	    if (n<=1)
		{
	       cout<<endl<<"节点数不符，需重新输入";     //确保节点数大于1
		}
	  }
	  double * x  = new double[n];
	  double * y  = new double[n];
      double * h = new double[n];
      
	  double * dy = new double[n];

	  cout<<"结点"<<ends<<ends<<ends<<"函数值"<<endl;
	  for(int m=0;m<=n;m++)
	  {
		x[m]=m*6.0/n;                                //插值节点x值
		y[m] = fun_func(x[m]);                 		 //插值节点函数值y
		cout<<x[m]<<ends<<ends<<ends<<ends<<ends<<y[m]<<endl;
	  }
	   while(v>6||v<0)
	  {
	     cout<<"input the value of the variable:";     //获得差值点
	     cin>>v;
	     if(v>6||v<0)                                  //确定差值点属于区间[0，6]
		 {
            cout<<"差值点不正确，请重新输入"<<endl;     
		 }
	   }                                        
	  

	  int i=judge(x,v,n);                            
	  double Result=Result_TriSpline(y, x,h,dy,v,i,n);
	  /*
	  delete []x;
	  delete []y;
	  delete []h;
	  delete []dy;
	  */
	  cout<<endl<<"三次样条插值所得函数值为:"<<Result<<endl;

	  double True_functionvalue = fun_func(v);

	  cout<<endl<<" 实  际  函  数  值  为 :"<<True_functionvalue<<endl;

	  cout<<endl<<"是否继续计算?[ture(1),false(0)]:"<<endl;

	  cin>>ToF_goon;

	}
//	return 0;
}