text.cpp

这是样条插值的源代码

#include<stdio.h>
#define MAX_N 20
typedef struct tagPOINT
{
   double x;
   double y;
}POINT;
void cha(POINT points[],double M[],double x,int n)
{
   	double S,q,p;
	int i,k;
	double h[MAX_N+1],b[MAX_N+1],c[MAX_N+1],d[MAX_N+1];
	double u[MAX_N+1],v[MAX_N+1],y[MAX_N+1];
	
	
	
	//计算M关系式中个参数的值
	h[0]=points[1].x-points[0].x;
	for(i=1;i<n;i++)
	{
	   h[i]=points[i+1].x-points[i].x;
	   b[i]=h[i]/(h[i]+h[i-1]);
	   c[i]=1-b[i];
	   d[i]=6*((points[i+1].y-points[i].y)/h[i]-(points[i].y-points[i-1].y)/h[i-1])/(h[i]+h[i-1]);

	}
	//用追赶法计算Mi,i=1,....n-1
	d[1]-=c[1]*M[0];
	d[n-1]-=b[n-1]*M[n];
	b[n-1]=0;c[1]=0;v[0]=0;
	for(i=1;i<n;i++)
	{
	  u[i]=2-c[i]*v[i-1];
	  v[i]=b[i]/u[i];
	  y[i]=(d[i]-c[i]*y[i-1])/u[i];

	}
	for(i=1;i<n;i++)
	{
	  M[n-i]=y[n-i]-v[n-i]*M[n-i+1];

	}
	//计算三次样条插值函数在X处的值
	k=0;
	while(x>=points[k].x)k++;
	k=k-1;
	p=points[k+1].x-x;
	q=x-points[k].x;
    S=(p*p*p*M[k]+q*q*q*M[k+1])/(6*h[k])+(p*points[k].y+q*points[k+1].y)/h[k]-h[k]*(p*M[k]+q*M[k+1])/6;
	printf("S(%f)=%f\n",x,S);


}
int main()
{
	int n;
	int i;
	POINT points[MAX_N+1];
	double M[MAX_N+1];
	
	double x;
	printf("\n Input n value:"); //输入插值点的数目
    scanf("%d",&n);
	if(n>MAX_N)
	{
	
		printf("The input n is larger the MAX_N,please redefine the MAX_N.\n");
		return 1;
	
	}

	if(n<=0)
	{
	  printf("Please input a number between 1 and %d.\n",MAX_N);
	  return 1;
	}
    //输入插值点（x_i,y_i),M0值和Mn值
    printf("Now input the(x_i,y_i),i=0,...,%d\n",n);
	for(i=0;i<=n;i++)
		scanf("%lf%lf",&points[i].x,&points[i].y);
	printf("Now input the M[0] value:");
	scanf("%lf",&M[0]);
	printf("Now input the M[n] value:");
	scanf("%lf",&M[n]);
	printf("Now input the x value:");//
	scanf("%lf",&x);
	if(x>points[n].x||x<points[0].x)
	{
		printf("Please input a number between %f and %f.\n",points[0].x,points[n].x);
		return 1;
	
	}


    cha(points,M,x,n);

	return 0;

}