scyt.cpp

三次样条插值函数

#include<iostream>
#include<iomanip>
using namespace std;

double zhuigan(double c[],double a[],double f[],double m0[],int n)
{	int i;
	double B[11],Y[12],X[12],b[12];
	for(i=0;i<12;i++) b[i]=2;
	B[0]=c[0]/b[0];
	for(i=1;i<11;i++) B[i]=c[i]/(b[i]-a[i-1]*B[i-1]);
	Y[0]=f[0]/b[0];
	for(i=1;i<12;i++) Y[i]=(f[i]-a[i-1]*Y[i-1])/(b[i]-a[i-1]*B[i-1]);
	X[n-1]=Y[n-1]; 
	for(i=n-2;i>=0;i--) X[i]=Y[i]-B[i]*X[i+1];
	for(i=0;i<12;i++) {cout<<"M"<<(i+1)<<" = "<<setiosflags(ios::scientific)<<X[i]<<"  "<<endl;m0[i]=X[i];}
	cout<<endl;
	return 0;
}

double yangtiao(double a,double b,double x[],double y[],double m[],int n)
{	int i;
	double h[11],shang[11],xia[11],d[12];
	for(i=0;i<n-1;i++)
		h[i]=x[i+1]-x[i];
	shang[0]=1;xia[10]=1;
	for(i=0;i<n-1;i++) shang[i]=h[i]/(h[i-1]+h[i]);
	for(i=0;i<n-1;i++) xia[i]=h[i]/(h[i]+h[i+1]);
	d[0]=(6/h[0])*(((y[1]-y[0])/h[0])-a);
    d[n-1]=(6/h[n-2])*(b-((y[n-1]-y[n-2])/h[n-2]));	
	for(i=1;i<n-1;i++) d[i]=(6/(h[i-1]+h[i]))*(((y[i+1]-y[i])/h[i])-((y[i]-y[i-1])/h[i-1]));
	zhuigan(shang,xia,d,m,n);
	return 0;
}

double chazhi(double x[],double y[],double m[],double t[],int n,int l)
{
	double h[11],z0[8],z1[8],z2[8];
	int i,j,k;
	for(i=0;i<n-1;i++)
		h[i]=x[i+1]-x[i];
	cout<<"节点"<<"              "<<"函数值"<<"              "<<"一阶导数"<<"              "
		<<"二阶导数"<<endl;
	for(i=0;i<l;i++)
	{	
		k=0;
		for(j=1;j<n;j++) {if(t[i]>x[j]) k++;}
		z0[i]=(m[k]/(6*h[k]))*(x[k+1]-t[i])*(x[k+1]-t[i])*(x[k+1]-t[i])
			+(m[k+1]/(6*h[k]))*(t[i]-x[k])*(t[i]-x[k])*(t[i]-x[k])
			+(y[k]-(m[k]*h[k]*h[k]/6))*((x[k+1]-t[i])/h[k])
			+(y[k+1]-(m[k+1]*h[k]*h[k]/6))*((t[i]-x[k])/h[k]);
		z1[i]=-(m[k]/(2*h[k]))*(x[k+1]-t[i])*(x[k+1]-t[i])
			 +(m[k+1]/(2*h[k]))*(t[i]-x[k])*(t[i]-x[k])
			 +(y[k+1]-y[k])/h[k]-((m[k+1]-m[k])/6)*(h[k]);
		z2[i]=(m[k]*(x[k+1]-t[i]))/h[k]+(m[k+1]*(t[i]-x[k]))/h[k];
		cout<<setiosflags(ios::scientific)<<t[i]<<"    "<<z0[i]<<"      "<<z1[i]<<"      "<<z2[i]<<"      "<<endl;
	}
	cout<<endl;
	return 0;
}

int main()
{
	double x0[12]={0.52,8.0,17.95,28.65,50.65,104.6,156.6,260.7,364.4,468.0,507.0,520.0};
	double y0[12]={5.28794,13.84,20.2,24.9,31.1,36.5,36.6,31.0,20.9,7.8,1.5,0.2};
    double t0[8]={4.0,14.0,30.0,60.0,130.0,230.0,450.0,515.0};
	double f0,fn;
	double M[12];
	f0=1.86548;fn=-0.046115;
	yangtiao(f0,fn,x0,y0,M,12);
	cout<<endl;
	chazhi(x0,y0,M,t0,12,8);
	return 0;
}