/*   三次樣條插值計算算法   */   
  #include   "math.h"   
  #include   "stdio.h"   
  #include   "stdlib.h"   
    
    
  /*   
  N：已知節(jié)點數(shù)N+1   
  R：欲求插值點數(shù)R+1   
  x，y為給定函數(shù)f(x)的節(jié)點值{x(i)}   (x(i)<x(i+1))   ,以及相應(yīng)的函數(shù)值{f(i)}     0<=i<=N   
  P0=f(x0)的二階導(dǎo)數(shù)；Pn=f(xn)的二階導(dǎo)數(shù)   
  u:存插值點{u(i)}       0<=i<=R   
  求得的結(jié)果s(ui)放入s[R+1]       0<=i<=R   
  返回0表示成功，1表示失敗   
  */   
  int   SPL(int   N,int   R,double   x[],double   y[],double   P0,double   Pn,double   u[],double   s[])   
  {   
  /*聲明局部變量*/   
  double   *h;   /*存放步長:{hi}       0<=i<=N-1     */   
  double   *a;   /*存放系數(shù)矩陣{ai}       1<=i<=N   ；   分量0沒有利用     */   
  double   *c;   /*先存放系數(shù)矩陣{ci}     后存放{Bi}     0<=i<=N-1     */   
  double   *g;   /*先存放方程組右端項{gi}     后存放求解中間結(jié)果{yi}     0<=i<=N     */   
  double   *af;   /*存放系數(shù)矩陣{a(f)i}       1<=i<=N   ；     */   
  double   *ba;   /*存放中間結(jié)果   0<=i<=N-1*/   
  double   *m;   /*存放方程組的解{m(i)}       0<=i<=N   ；     */   
    
  int   i,k;   
  double   p1,p2,p3,p4;       
    
  /*分配空間*/   
  if(!(h=(double*)malloc(N*sizeof(double))))   exit(1);   
  if(!(a=(double*)malloc((N+1)*sizeof(double))))   exit(1);   
  if(!(c=(double*)malloc(N*sizeof(double))))   exit(1);   
  if(!(g=(double*)malloc((N+1)*sizeof(double))))   exit(1);   
  if(!(af=(double*)malloc((N+1)*sizeof(double))))   exit(1);   
  if(!(ba=(double*)malloc((N)*sizeof(double))))   exit(1);   
  if(!(m=(double*)malloc((N+1)*sizeof(double))))   exit(1);   
    
  /*第一步：計算方程組的系數(shù)*/   
  for(k=0;k<N;k++)   
  h[k]=x[k+1]-x[k];   
  for(k=1;k<N;k++)   
  a[k]=h[k]/(h[k]+h[k-1]);   
  for(k=1;k<N;k++)   
  c[k]=1-a[k];   
  for(k=1;k<N;k++)   
  g[k]=3*(c[k]*(y[k+1]-y[k])/h[k]+a[k]*(y[k]-y[k-1])/h[k-1]);   
  c[0]=a[N]=1;   
  g[0]=3*(y[1]-y[0])/h[0]-P0*h[0]/2;   
  g[N]=3*(y[N]-y[N-1])/h[N-1]+Pn*h[N-1]/2;   
    
  /*第二步：用追趕法解方程組求{m(i)}     */   
  ba[0]=c[0]/2;   
  g[0]=g[0]/2;   
  for(i=1;i<N;i++)         
  {   
  af[i]=2-a[i]*ba[i-1];   
  g[i]=(g[i]-a[i]*g[i-1])/af[i];   
  ba[i]=c[i]/af[i];   
  }   
  af[N]=2-a[N]*ba[N-1];   
  g[N]=(g[N]-a[N]*g[N-1])/af[N];   
    
  m[N]=g[N];                   /*P110   公式：6.32*/   
  for(i=N-1;i>=0;i--)   
  m[i]=g[i]-ba[i]*m[i+1];   
    
  /*第三步：求值*/   
  for(i=0;i<=R;i++)   
  {   
  /*判斷u(i)屬于哪一個子區(qū)間，即確定k   */   
  if(u[i]<x[0]   ||   u[i]>x[N])       
  {   
  /*釋放空間*/   
  free(h);   
  free(a);   
  free(c);   
  free(g);   
  free(af);   
  free(ba);   
  free(m);   
  return   1;   
  }   
  k=0;   
  while(u[i]>x[k+1])   
  k++;   
    
    
  p1=(h[k]+2*(u[i]-x[k])*pow((u[i]-x[k+1]),2)*y[k])/pow(h[k],3);   
  p2=(h[k]-2*(u[i]-x[k+1])*pow((u[i]-x[k]),2)*y[k+1])/pow(h[k],3);   
  p3=(u[i]-x[k])*pow((u[i]-x[k+1]),2)*m[k]/pow(h[k],2);   
  p4=(u[i]-x[k+1])*pow((u[i]-x[k]),2)*m[k+1]/pow(h[k],2);   
  s[i]=p1+p2+p3+p4;   
  }   
    
  /*釋放空間*/   
  free(h);   
  free(a);   
  free(c);   
  free(g);   
  free(af);   
  free(ba);   
  free(m);   
    
  return   0;   
  }   
    
    
  void   main()   
  {   
  int   N,R;   
  double   *x,*y,*u,*s;   
  double   P0,Pn;   
  int   i;   
    
  /*驗證算法：*/   
  N=7;   
  R=6;   
    
  /*分配空間*/   
  if(!(x=(double*)malloc((N+1)*sizeof(double))))   
  {   
  printf("malloc   error!\n");   
  exit(1);   
  }   
  if(!(y=(double*)malloc((N+1)*sizeof(double))))   
  {   
  printf("malloc   error!\n");   
  exit(1);   
  }   
  if(!(u=(double*)malloc((R+1)*sizeof(double))))   
  {   
  printf("malloc   error!\n");   
  exit(1);   
  }   
  if(!(s=(double*)malloc((R+1)*sizeof(double))))   
  {   
  printf("malloc   error!\n");   
  exit(1);   
  }   
    
  x[0]=0.5;x[1]=0.7;x[2]=0.9;x[3]=1.1;x[4]=1.3;x[5]=1.5;x[6]=1.7;x[7]=1.9;   
  y[0]=0.4794;y[1]=0.6442;y[2]=0.7833;y[3]=0.8912;y[4]=0.9636;y[5]=0.9975;y[6]=0.9917;y[7]=0.9463;   
  u[0]=0.6;u[1]=0.8;u[2]=1.0;u[3]=1.2;u[4]=1.4;u[5]=1.6;u[6]=1.8;   
    
  P0=-0.4794;   
  Pn=-0.9463;   
    
  if(!SPL(   N,   R,   x,   y,   P0,   Pn,   u,   s))   
  {   
  /*打印結(jié)果*/   
  printf("\nx=   ");   
  for(i=0;i<=N;i++)   
  printf("%8.1f",x[i]);   
    
  printf("\ny=   ");   
  for(i=0;i<=N;i++)   
  printf("%8.4f",y[i]);   
    
  printf("\n\nu=       ");   
  for(i=0;i<=R;i++)   
  printf("%9.2f",u[i]);   
  printf("\ns=       ");   
  for(i=0;i<=R;i++)   
  printf("%9.5f",s[i]);   
  printf("\nsin=   ");   
  for(i=0;i<=R;i++)   
  printf("%9.5f",sin(u[i]));   
  }   
    
  /*釋放空間*/   
  free(x);   
  free(y);   
  free(u);   
  free(s);   
  }
Top