insert.cpp

三次样条插值

#include"iostream.h"
#include"math.h"
#include"stdlib.h"

amspl(double x[], double y[], int n, double dy[], double ddy[], double t[], int m, double z[], double dz[], double ddz[]);

void main()
{
	const int size=5;
	int n=size;
	const int m=200;
	double s;
	double x[size]={0,50,100,150,200};
		//,150,100,50};
	double y[size]={50,20,0,20,50};
		//,80,100,80};
	double dy[size],ddy[size],t[m],z[m],dz[m],ddz[m];
	for(int i=0;i<m;i++)
	{
		t[i]=i;
	}

/*	for( i=m/2;i<m;i++)
	{
		t[i]=m-i;
	}
*/
	dy[0]=-1.0;
	dy[size-1]=1.0;
	s=amspl(x,y,n,dy,ddy,t,m,z,dz,ddz);
	
	for( i=0;i<m;i++)
	{
		cout<<"z["<<i<<"]="<<z[i]<<",";
		if(i%5==0)
		{
			cout<<endl;
		}
	}
		
	
}
amspl(double x[], double y[], int n, double dy[], double ddy[], double t[], int m, double z[], double dz[], double ddz[])
{
	int i,j;
	double h0,h1,alpha,deta,g;
	double *s;
	s=(double *) malloc(n*sizeof(double));
	s[0]=dy[0];
	dy[0]=0.0;
	h0=x[1]-x[0];
	for(j=1;j<=n-2;j++)
	{
		h1=x[j+1]-x[j];
		alpha=h0/(h0+h1);
		deta=(1.0-alpha)*(y[j]-y[j-1])/h0;
		deta=3.0*(deta+alpha*(y[j+1]-y[j])/h1);
		dy[j]=-alpha/(2.0+(1.0-alpha)*dy[j-1]);
		s[j]=(deta-(1.0-alpha)*s[j-1]);
		s[j]=s[j]/(2.0+(1.0-alpha)*dy[j-1]);
		h0=h1;
	}
	for(j=n-2;j>=0;j--)
	 dy[j]=dy[j]*dy[j+1]+s[j];
	for(j=0;j<=n-2;j++)
	s[j]=x[j+1]-x[j];
	for(j=0;j<=n-2;j++)
	 {
		h1=s[j]*s[j];
		ddy[j]=6.0*(y[j+1]-y[j])/h1-2.0*(2.0*dy[j]+dy[j+1])/s[j];
	}
	h1=s[n-2]*s[n-2];
	ddy[n-1]=6.0*(y[n-2]-y[n-1])/h1+2.0*(2.0*dy[n-1]+dy[n-2])/s[n-2];
	g=0.0;
	for(i=0;i<=n-2;i++)
	{
		h1=0.5*s[i]*(y[i]+y[i+1]);
		h1=h1-s[i]*s[i]*s[i]*(ddy[i]+ddy[i+1])/24.0;
		g=g+h1;
	}
	for(j=0;j<=m-1;j++)
	{
		if (t[j]>=x[n-1])
			i=n-2;
		else
		{
			i=0;
			while(t[j]>x[i+1])
				i=i+1;
		}
		h1=(x[i+1]-t[j])/s[i];
		h0=h1*h1;
		z[j]=(3.0*h0-2.0*h0*h1)*y[i];
		z[j]=z[j]+s[i]*(h0-h0*h1)*dy[i];
		dz[j]=6.0*(h0-h1)*y[i]/s[i];
		dz[j]=dz[j]+(3.0*h0-2.0*h1)*dy[i];
		ddz[j]=(6.0-12.0*h1)*y[i]/(s[i]*s[i]);
		ddz[j]=ddz[j]+(2.0-6.0*h1)*dy[i]/s[i];
		h1=(t[j]-x[i])/s[i];
		h0=h1*h1;
		z[j]=z[j]+(3.0*h0-2.0*h0*h1)*y[i+1];
		z[j]=z[j]-s[i]*(h0-h0*h1)*dy[i+1];
		dz[j]=dz[j]-6.0*(h0-h1)*y[i+1]/s[i];
		dz[j]=dz[j]+(3.0*h0-2.0*h1)*dy[i+1];
		ddz[j]=ddz[j]+(6.0-12.0*h1)*y[i+1]/(s[i]*s[i]);
		ddz[j]=ddz[j]-(2.0-6.0*h1)*dy[i+1]/s[i];
		
	}
delete(s);

return(g);
}