自然样条函数.cpp

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

#include<stdio.h>
#define MAX_N 20           // 定义（x_i,y_i)的最大维数
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;
 }
    //输入插值点（x_i,y_i) 
 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);
 M[0]=0,M[n]=0;
 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;
	 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
 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];
   //计算插值函数在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(%lf)=%f\n",x,S);
 return 0;
}