注册 登录  
 加关注
查看详情
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

回首望星辰

See you in the next world

 
 
 

日志

 
 

三次样条插值c代码  

2008-11-26 14:12:37|  分类: 图形图像开发 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

///////////////////////////////////////////////////

// purpose:给定,值的三次样条插值多项式  //

///////////////////////////////////////////////////
# 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),M0值和Mn值
 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("%lf %lf\n",points[i].x,points[i].y);
 }
 printf ("Now input the M[0] value:");
 scanf ("%lf",&M[0]);
 printf ("\nNow input the M[n] value:");
 scanf ("%lf",&M[n]);
 printf ("\nNow input the x value:");    
 //输入计算三次样条插值函数的x值
 scanf ("%lf",&x);
 if (x>points[n-1]. x || x<points[0].x)
 {
  printf ("Please input a number between %f and %f. \n",points[0].x,points[n].x);
  return 1;
 }
 // 计算M关系式中各参数的值
 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,i=1,…,n-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=1;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]+q* M[k+1])/6;
 printf ("S(%f ) = %f\n",x,S);
 // 输出
 getchar ( );
 return 0;
}

  评论这张
 
阅读(1247)| 评论(0)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2018