三次样条插值函数(C++版)

  1. double spl2(double x[],                                                                /*x坐标序列*/
  2.                         double y[],                                                                /*y坐标序列*/
  3.                         int n,                                                                        /*输入数据个数*/
  4.                         double ddy1, double ddyn,                                /*第一点和最末点二阶导数*/
  5.                         double t[],                                                                /*插值点的x坐标序列*/
  6.                         int m,                                                                        /*插值点个数*/
  7.                         double z[]                                                                /*差值点的y坐标序列*/
  8. )        
  9. {
  10.         int i,j;
  11.         double h0,h1,alpha,beta,*s,*dy;
  12.         s= new double[n];
  13.         dy = new double[n];
  14.         dy[0]=-0.5;
  15.         h0=x[1]-x[0];
  16.         s[0]=3.0*(y[1]-y[0])/(2.0*h0)-ddy1*h0/4.0;
  17.         for (j=1;j<=n-2;j++)
  18.         {
  19.                 h1=x[j+1]-x[j];
  20.                 alpha=h0/(h0+h1);
  21.                 beta=(1.0-alpha)*(y[j]-y[j-1])/h0;
  22.                 beta=3.0*(beta+alpha*(y[j+1]-y[j])/h1);
  23.                 dy[j]=-alpha/(2.0+(1.0-alpha)*dy[j-1]);
  24.                 s[j]=(beta-(1.0-alpha)*s[j-1]);
  25.                 s[j]=s[j]/(2.0+(1.0-alpha)*dy[j-1]);
  26.                 h0=h1;
  27.         }
  28.         dy[n-1]=(3.0*(y[n-1]-y[n-2])/h1+ddyn*h1/2.0-s[n-2])/(2.0+dy[n-2]);
  29.         for (j=n-2;j>=0;j--)        dy[j]=dy[j]*dy[j+1]+s[j];
  30.         for (j=0;j<=n-2;j++)        s[j]=x[j+1]-x[j];
  31.         for (j=0;j<=m-1;j++)
  32.         {
  33.                 if (t[j]>=x[n-1]) i=n-2;
  34.                 else
  35.                 {
  36.                         i=0;
  37.                         while (t[j]>x[i+1]) i=i+1;
  38.                 }
  39.                 h1=(x[i+1]-t[j])/s[i];
  40.                 h0=h1*h1;
  41.                 z[j]=(3.0*h0-2.0*h0*h1)*y[i];
  42.                 z[j]=z[j]+s[i]*(h0-h0*h1)*dy[i];
  43.                 h1=(t[j]-x[i])/s[i];
  44.                 h0=h1*h1;
  45.                 z[j]=z[j]+(3.0*h0-2.0*h0*h1)*y[i+1];
  46.                 z[j]=z[j]-s[i]*(h0-h0*h1)*dy[i+1];
  47.         }
  48.         delete[] s;
  49.         delete[] dy;
  50.         return 0.0;
  51. }
  • 4
    点赞
  • 23
    收藏
    觉得还不错? 一键收藏
  • 8
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 8
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值