spline.txt

三次样条插值的源代码

    
  /*   三次样条插值计算算法   */   
  #include   "math.h"   
  #include   "stdio.h"   
  #include   "stdlib.h"   
    
    
  /*   
  N：已知节点数N+1   
  R：欲求插值点数R+1   
  x，y为给定函数f(x)的节点值{x(i)}   (x(i)<x(i+1))   ,以及相应的函数值{f(i)}     0<=i<=N   
  P0=f(x0)的二阶导数；Pn=f(xn)的二阶导数   
  u:存插值点{u(i)}       0<=i<=R   
  求得的结果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;   /*存放系数矩阵{ai}       1<=i<=N   ；   分量0没有利用     */   
  double   *c;   /*先存放系数矩阵{ci}     后存放{Bi}     0<=i<=N-1     */   
  double   *g;   /*先存放方程组右端项{gi}     后存放求解中间结果{yi}     0<=i<=N     */   
  double   *af;   /*存放系数矩阵{a(f)i}       1<=i<=N   ；     */   
  double   *ba;   /*存放中间结果   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);   
    
  /*第一步：计算方程组的系数*/   
  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)属于哪一个子区间，即确定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))   
  {   
  /*打印结果*/   
  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