三次样条_070460.cpp

对于给定（n+1)个节点上的函数值

#include<iostream.h>
#include<math.h>
#include<fstream.h>

//输入节点位置及函数值
void init_xy_1(int* n,double x[100],double y[100],double* Do,double* Dn)
{
	double d0,dn;
	int m;
	cout<<"请输入插值节点个数N:";
	cin>>m;
	*n=m-1;
    cout<<"请输入这"<<m<<"个节点的位置Xi(i=0,1,..."<<m-1<<"):"<<endl;
	for(int i=0;i<=m-1;i++)
	cin>>x[i];
	cout<<"请输入这"<<m<<"个节点的函数值Yi(i=0,1,..."<<m-1<<"):"<<endl;
	for(int j=0;j<=m-1;j++)
	cin>>y[j];
	cout<<"请输入一型插值条件"<<endl;
	cout<<"第0个节点处的一阶导数:";
	cin>>d0;
	*Do=d0;
    cout<<"第"<<m-1<<"个节点处的一阶导数:";
	cin>>dn;
	*Dn=dn;
	cout<<endl;

}



//计算正规方程组所需参数
void u_v_d(int n,double x[100],double y[100],double u[100],double v[100],double d[100],double Do,double Dn)
{
	u[0]=0;
	v[0]=0;
	for(int i=1;i<=n-1;i++)
	{
		u[i]=(double)(x[i]-x[i-1])/(x[i+1]-x[i-1]);
		v[i]=(double)(x[i+1]-x[i])/(x[i+1]-x[i-1]);
		d[i]=(double)6.0*((y[i+1]-y[i])/(x[i+1]-x[i])-(y[i]-y[i-1])/(x[i]-x[i-1]))/(x[i+1]-x[i-1]);
	};
	d[0]=(double)6*((y[1]-y[0])/(x[1]-x[0])-Do)/(x[1]-x[0]);
	d[n]=(double)6*(Dn-(y[n]-y[n-1])/(x[n]-x[n-1]))/(x[n]-x[n-1]);
}

//初始化正规方程组
void init_mat(int n,double u[100],double v[100],double d[100],double matrix[][100])
{
	for(int i=0;i<=n;i++)
		for(int j=0;j<=n;j++)
			matrix[i][j]=0;
	for(int l=0;l<=n;l++)
		matrix[l][l]=2;
	for(int k=1;k<=n-1;k++)
	{
		matrix[k][k-1]=u[k];
		matrix[k][k+1]=v[k];
	}
	matrix[0][1]=1;
	matrix[n][n-1]=1;
}

//追赶法解正规方程组，求出M值
void caculate(int n,double a[][100],double d[100],double m[100])
{
	int i,j;
	double li;
	for(i=1;i<=n;i++)
    {
		li=(double)a[i][i-1]/a[i-1][i-1];
	  for(j=i-1;j<=n;j++)
	  {
		  a[i][j]=a[i][j]-li*a[i-1][j];
	  }
      d[i]=d[i]-d[i-1]*li;
	}
	m[n]=(double)d[n]/a[n][n];
	for(i=n-1;i>=0;i--)
	{
		m[i]=(double)(d[i]-a[i][i+1]*m[i+1])/a[i][i];
	}
	  
}

void selectprint(double y)
{
	if ((y>0)&&(y!=1)) cout<<"+"<<y;
      else if (y==1) cout<<"+";
          else if (y<0) cout<<y;
              
}

//输出S(x)的表达式
void print(int n,double x[100],double y[100],double m[100],double e0[100],double e1[100],double e2[100],double e3[100])
{
	
	for(int i=0;i<=n;i++)
	{
		e0[i]=y[i];
		e1[i]=(double)((double)(y[i+1]-y[i])/(x[i+1]-x[i])-(1.0/3*m[i]+1.0/6*m[i+1])*(x[i+1]-x[i]));
		e2[i]=(double)1.0/2*m[i];
		e3[i]=(double)(1.0/(6*(x[i+1]-x[i])))*(m[i+1]-m[i]);
		if((i>=0)&&(i<=n-1))
		{
		cout<<"当X值在区间["<<x[i]<<","<<x[i+1]<<"]时"<<endl;
		cout<<"S(x)=";
	
		if(e0[i])
		cout<<e0[i];
		if(e1[i])
		{
			selectprint(e1[i]);
		    cout<<"(x-"<<x[i]<<")";
		
		}
		if(e2[i])
		{
			selectprint(e2[i]);
            cout<<"(x-"<<x[i]<<")"<<"(x-"<<x[i]<<")";
		
		}
		if(e3[i])
		{
			selectprint(e3[i]);
            cout<<"(x-"<<x[i]<<")"<<"(x-"<<x[i]<<")"<<"(x-"<<x[i]<<")"<<'\n';
		}
		cout<<endl;
		}
	}
}

//计算本题要求打印的节点函数值
void prints_1(int n,double x[100],double e0[100],double e1[100],double e2[100],double e3[100])
{
	double s[100];
	cout<<"打印S(i+0.5)的值"<<endl;
	for(int i=0;i<=n-1;i++)
	{
        s[i]=(double)e0[i]+e1[i]*(i+0.5-x[i])+e2[i]*(i+0.5-x[i])*(i+0.5-x[i])+e3[i]*(i+0.5-x[i])*(i+0.5-x[i])*(i+0.5-x[i]);
		cout<<"S("<<i+0.5<<")="<<s[i]<<endl;
	}
}

//输出任意点的函数值

void prints_2(int n,double x[100],double e0[100],double e1[100],double e2[100],double e3[100])
{
	double en,sn;
	int ki;
	cout<<"输入任意一个节点的位置:"<<"(所取范围为:["<<x[0]<<","<<x[n]<<"])";
	cin>>en;
	if(en>=x[0]&&en<=x[n])
	{
	do{
		for(int j=0;j<=n;j++)
		{
		if(x[j]<=en)
			ki=j;
		}

        sn=(double)e0[ki]+e1[ki]*(en-x[ki])+e2[ki]*(en-x[ki])*(en-x[ki])+e3[ki]*(en-x[ki])*(en-x[ki])*(en-x[ki]);
		cout<<"S("<<en<<")="<<sn<<endl;
        cout<<"输入下一个所求点的位置:"<<"(所取范围为:["<<x[0]<<","<<x[n]<<"])";
		cin>>en;
	}
    while(en>=x[0]&&en<=x[n]);
	}
	else cout<<"输入节点不在取值范围内!"<<endl;
}


     
void main()
{
  	double x[100],y[100],u[100],v[100],d[100],m[100];
    double e0[100],e1[100],e2[100],e3[100];//表达式系数
	double matrix[100][100];
	double Do,Dn;
	int n,an;

	ofstream outfile;
	ifstream infile;
    outfile.open("print.txt",ios::trunc);
	infile.open("jiedian.txt",ios::in);
	
    cout<<"若计算书中所求函数_请选择1:"<<endl;
	cout<<"求其它任意3次样条插值函数_请选择2:"<<endl;

	cin>>an;

	if(an==1)
	{	
		for(int j=0;j<=10;j++)
		{
		   infile>>x[j];
		   infile>>y[j];
		}
    	infile>>Do;
	    infile>>Dn;
		n=10;

	    u_v_d(n,x,y,u,v,d,Do,Dn);
	
	    init_mat(n,u,v,d,matrix);
	
        caculate(n,matrix,d,m);

	    print(n,x,y,m,e0,e1,e2,e3);

        prints_1(n,x,e0,e1,e2,e3);


		for(int i=0;i<=n-1;i++)
		{
          outfile<<"当X值在区间["<<x[i]<<","<<x[i+1]<<"]时"<<endl;
	      outfile<<"S(x)=";
		
		  if(e0[i])
	      outfile<<e0[i];
		  if(e1[i])
		  {
			if ((e1[i]>0)&&(e1[i]!=1)) outfile<<"+"<<e1[i];
              else if (e1[i]==1) outfile<<"+";
                else if (e1[i]<0) outfile<<e1[i];
		    
			outfile<<"(x-"<<x[i]<<")";
		  }
		  if(e2[i])
		  {
			 if ((e2[i]>0)&&(e2[i]!=1)) outfile<<"+"<<e2[i];
              else if (e2[i]==1) outfile<<"+";
                else if (e2[i]<0) outfile<<e2[i];
		    
				outfile<<"(x-"<<x[i]<<")"<<"(x-"<<x[i]<<")";
		  }
		  if(e3[i])
		  {
			if ((e3[i]>0)&&(e3[i]!=1)) outfile<<"+"<<e3[i];
              else if (e3[i]==1) outfile<<"+";
                else if (e3[i]<0) outfile<<e3[i];
		    
            
			outfile<<"(x-"<<x[i]<<")"<<"(x-"<<x[i]<<")"<<"(x-"<<x[i]<<")"<<'\n';
		  }
		
		  outfile<<endl;
		}
		outfile<<"打印S(i+0.5)的值"<<endl;
	    for(int l=0;l<=n-1;l++)
		{
        outfile<<"S("<<l+0.5<<")="<<e0[l]+e1[l]*(l+0.5-x[l])+e2[l]*(l+0.5-x[l])*(l+0.5-x[l])+e3[l]*(l+0.5-x[l])*(l+0.5-x[l])*(l+0.5-x[l])<<endl;
		}
		prints_2(n,x,e0,e1,e2,e3);
	}
	
	if(an==2)
	{
	    init_xy_1(&n,x,y,&Do,&Dn);

	    u_v_d(n,x,y,u,v,d,Do,Dn);
	
	    init_mat(n,u,v,d,matrix);
	
        caculate(n,matrix,d,m);

    	print(n,x,y,m,e0,e1,e2,e3);

        prints_2(n,x,e0,e1,e2,e3);
	}


infile.close();
outfile.close();
}