m型自然样条插值.cpp

利用C语言编写的一些常见的差值程序

#include<stdio.h>
#define MAX_N 20
typedef struct tagPOINT
{double x;
 double y;
}POINT;
int main()
{int n;
 int i,k;
 POINT points[MAX_N+1];
 double h[MAX_N+1],b[MAX_N+1],c[MAX_N+1],d[MAX_N+1],m[MAX_N+1];
 double u[MAX_N+1],v[MAX_N+1],y[MAX_N+1];
 double x,p,q,S;
 printf("\nInput n value:");
 scanf("%d",&n);
 if(n>MAX_N)
 {printf("The input n is larger than 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;
 }
 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] vale:");
     scanf("%f",&m[0]);
 printf("Now input the m[n] vale:");
     scanf("%f",&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 %lf and %lf.\n",points[0].x,points[n].x);
  return 1;
 }
 h[0]=points[1].x-points[0].x;
 for(i=1;i<n;i++)
 {
	 h[i]=points[i+1].x-points[i].x;
	 c[i]=h[i]/(h[i]+h[i-1]);
	 b[i]=1-c[i];
	 d[i]=3*(c[i]*(points[i+1].y-points[i].y)/h[i]+b[i]*(points[i].y-points[i-1].y)/h[i-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=0;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];
 k=0;
 while(x>=points[k].x)
	 k++;
 k=k-1;
 p=points[k+1].x-x;
 q=x-points[k].x;
 S=(1+2*q/h[k])*p*p*points[k].y/(h[k]*h[k])+q*p*p*m[k]/(h[k]*h[k])+(1+2*p/h[k])*q*q*points[k+1].y/(h[k]*h[k])-p*q*q*m[k+1]/(h[k]*h[k]);
 printf("S(%lf)=%f\n",x,S);
 return 0;
}