报告二 2样条.cpp

样条插值函数

#include "stdlib.h"
#include "stdio.h"
#include "conio.h"
#include "math.h"
#define ERROR 0

int Locatexc(float *x,float xc,int n)//确定xc所在区间
{
 int i;
 for(i=0;i<=n;i++)
  {
   if((x[i]<=xc)&&xc<=x[i+1])
     return i;
  }
 return -1;
}

main()
{
 int n,i,loc;
 float *x,*y,*h,*a,*b,*A,*B,*m,xc,result;
 printf("请输入n(节点个数n+1):");
 scanf("%d",&n);
 x=(float*)malloc((n+1)*sizeof(float));
 y=(float*)malloc((n+1)*sizeof(float));
 for(i=0;i<=n;i++)
  {
   printf("请输入第%d个节点x值:\n",i);
   fflush(stdin);
   scanf("%f",&x[i]);
   printf("请输入第%d个节点y值:\n",i);
   scanf("%f",&y[i]);
  }
 printf("您输入的数据为:\ny = :");
 for(i=0;i<=n;i++)
  {
   printf("%f   ",y[i]);
   if(i==n)
   printf("\n");
  }
 printf("x = :");
 for(i=0;i<=n;i++)
  {
   printf("%f   ",x[i]);
   if(i==n)
   printf("\n");

  }

 loop:printf("请输入待求x的值:\n");
      fflush(stdin);
      scanf("%f",&xc);
      loc=Locatexc(x,xc,n);
      if(loc!=-1)printf("此x在x[%d]与x[%d]之间\n",loc,(loc+1));
       else
       {
        printf("此x不在这些区间中,请重新输入:\n");
        goto loop;
       }

  h=(float*)malloc(n*sizeof(float));
  for(i=0;i<=n-1;i++)
    h[i]=x[i+1]-x[i];

  a=(float*)malloc((n+1)*sizeof(float));
  a[0]=1;a[n]=0;
  for(i=1;i<=n-1;i++)
   a[i]=h[i-1]/(h[i-1]+h[i]);

  b=(float*)malloc((n+1)*sizeof(float));
  b[0]=3*(y[1]-y[0])/(h[0]); b[n]=3*(y[n]-y[n-1])/(h[n-1]);
  for(i=1;i<=n-1;i++)
   b[i]=3*((1-a[i])*(y[i]-y[i-1])/h[i-1]+a[i]*(y[i+1]-y[i])/h[i]);

  A=(float*)malloc(n*sizeof(float));
  A[0]=-a[0]/2;
  for(i=1;i<=n-1;i++)
   A[i]=-a[i]/(2+(1-a[i])*A[i-1]);

  B=(float*)malloc(n*sizeof(float));
  B[0]=b[0]/2;
  for(i=1;i<=n-1;i++)
   B[i]=(b[i]-(1-a[i])*B[i-1])/(2+(1-a[i])*A[i-1]);

  m=(float*)malloc((n+1)*sizeof(float));
  m[n]=(b[n]-(1-a[n])*B[n-1])/(2+(1-a[n])*A[n-1]);
  for(i=n-1;i>=0;i--)
   m[i]=A[i]*m[i+1]+B[i];

  result=(1+2*(xc-x[loc])/(x[loc+1]-x[loc]))*pow((xc-x[loc+1])/(x[loc]-x[loc+1]),2)*y[loc]+(1+2*(xc-x[loc+1])/(x[loc]-x[loc+1]))*pow((xc-x[loc])/(x[loc+1]-x[loc]),2)*y[loc+1]+(xc-x[loc])*pow((xc-x[loc+1])/(x[loc]-x[loc+1]),2)*m[loc]+(xc-x[loc+1])*pow((xc-x[loc])/(x[loc+1]-x[loc]),2)*m[loc+1];

  printf("函数在x点的近似值为%f\n",result);

  free(x);
  free(y);
  free(h);
  free(a);
  free(b);
  free(A);
  free(B);

 getch();
}