三次样条(第一类边值问题).txt

用于求解三次样条拟合！！！程序经过优化处理

#include<iostream.h>//三次样条插值(第一类边值问题)
double *Catch(double *a,double *b,double *c,double *f);
int n;
void main()
{	
	cout<<"请输入插值节点个数:"<<endl;
	cout<<"n=";
	cin>>n;
	double *X=new double [n];
	cout<<"请输入各节点的值:"<<endl;
	for(int i=0;i<n;i++)
	{
		cout<<"X["<<i<<"]=";
		cin>>X[i];
	}
	double *h=new double [n-1];
	for(int k=0;k<n-1;k++)
	h[k]=X[k+1]-X[k];
	double *a=new double [n-1];
	for(int j=0;j<n-2;j++)
		a[j]=h[j]/(h[j]+h[j+1]);
	a[n-2]=1;
	double *b=new double [n];
	for(int p=0;p<n;p++)
		b[p]=2;
	double *c=new double [n-1];
	c[0]=1;
	for(int q=1;q<n-1;q++)
		c[q]=1-a[q-1];
	double *f=new double [n];
	cout<<"请输入插值条件:"<<endl;
	for(int w=0;w<n;w++)
	{
		cout<<"f["<<X[w]<<"]=";
		cin>>f[w];
	}
	double *detaf=new double [n-1];
	for(int s=0;s<n-1;s++)
	detaf[s]=f[s+1]-f[s];
	double f1,f2;
	cout<<"请输入边界条件:"<<endl;
	cout<<"f'["<<X[0]<<"]=";
	cin>>f1;
	cout<<"f'["<<X[n-1]<<"]=";
	cin>>f2;
	double *d=new double [n];
	d[0]=6*(detaf[0]/h[0]-f1)/h[0];
	for(int m=1;m<n-1;m++)
		d[m]=6*(detaf[m]/h[m]-detaf[m-1]/h[m-1])/(h[m-1]+h[m]);
	d[n-1]=6*(f2-detaf[n-2]/h[n-2])/h[n-2];
	double *M=new double [n];
	M=Catch(a,b,c,d);
	cout<<"求得各节点处的二阶导数值为:"<<endl;
	for(int r=0;r<n;r++)
		cout<<"M["<<r<<"]="<<M[r]<<endl;
	cout<<"此插值问题的解为:"<<endl;
	for(int t=0;t<n-1;t++)
		cout<<"s(x)="<<M[t]/(6*h[t])<<"("<<X[t+1]<<'-'<<'x'<<")^3+"
			<<M[t+1]/(6*h[t])<<"(x-"<<X[t]<<")^3+"<<(f[t]-h[t]*h[t]*M[t]/6)/h[t]
			<<"("<<X[t+1]<<"-x)+"<<(f[t+1]-h[t]*h[t]*M[t+1]/6)/h[t]<<"(x-"<<X[t]<<"),"
			<<" "<<"x is in ["<<X[t]<<','<<X[t+1]<<']'<<endl;
	delete [] X;
	delete [] h;
	delete [] a;
	delete [] b;
	delete [] c;
	delete [] f;
	delete [] detaf;
	delete [] d;
	delete [] M;
}
double *Catch(double *a,double *b,double *c,double *f)
{
	double *s=new double [n];
	double *t=new double [n-1];
	s[0]=b[0];
	t[0]=c[0]/s[0];
	for(int i=1;i<n-1;i++)
	{
		s[i]=b[i]-a[i-1]*t[i-1];
		t[i]=c[i]/s[i];
	}
	s[n-1]=b[n-1]-a[n-2]*t[n-2];
	double *y=new double [n];
	y[0]=f[0]/s[0];
	for(int j=1;j<n;j++)
		y[j]=(f[j]-a[j-1]*y[j-1])/s[j];
	double *x=new double [n];
	x[n-1]=y[n-1];
	for(int k=n-2;k>=0;k--)
		x[k]=y[k]-t[k]*x[k+1];
	return x;
	delete [] s;
	delete [] t;
	delete [] y;
	delete [] x;
}