找回密码
 注册
Simdroid-非首页
查看: 50|回复: 0

[代码挑战区] 求助:三次样条插值问题

[复制链接]
发表于 2007-5-15 16:21:59 | 显示全部楼层 |阅读模式 来自 四川成都
使用三次样条插值对等间距的三个点进行插值(每两点间插100个点),y=(15.741419,14.706263,13.007044),从点的情况看应该是递减的一条曲线,生成出来的插值中15.741419这个值应该是第一个,但这个值却在100个点以后出现,这是为什么啊?三次样条插值代码如下:请大侠们帮忙看看!
double espl3(int * x,double * y,int n,int * t,int m, double * z)
{
    int i,j;
    double h0, y0, h1, y1, alpha, beta, u, g, *s, *dy, *ddy, *dz, *ddz;
    s=malloc(n*sizeof(double));
    dy=malloc(n*sizeof(double));
    ddy=malloc(n*sizeof(double));
    dz=malloc(m*sizeof(double));
    ddz=malloc(m*sizeof(double));
   
    h0=x[n-1]-x[n-2];
    y0=y[n-1]-y[n-2];
    dy[0]=0.0; ddy[0]=0.0; ddy[n-1]=0.0;
    s[0]=1.0; s[n-1]=1.0;
    for (j=1;j<=n-1;j++)
      { h1=h0; y1=y0;
        h0=x[j]-x[j-1];
        y0=y[j]-y[j-1];
        alpha=h1/(h1+h0);
        beta=3.0*((1.0-alpha)*y1/h1+alpha*y0/h0);
        if (j<n-1)
          { u=2.0+(1.0-alpha)*dy[j-1];
            dy[j]=-alpha/u;
            s[j]=(alpha-1.0)*s[j-1]/u;
            ddy[j]=(beta-(1.0-alpha)*ddy[j-1])/u;
          }
      }
    for (j=n-2;j>=1;j--)
      { s[j]=dy[j]*s[j+1]+s[j];
        ddy[j]=dy[j]*ddy[j+1]+ddy[j];
      }
    dy[n-2]=(beta-alpha*ddy[1]-(1.0-alpha)*ddy[n-2])/
            (alpha*s[1]+(1.0-alpha)*s[n-2]+2.0);
    for (j=2;j<=n-1;j++)
        dy[j-2]=s[j-1]*dy[n-2]+ddy[j-1];
    dy[n-1]=dy[0];
    for (j=0;j<=n-2;j++) s[j]=x[j+1]-x[j];
    for (j=0;j<=n-2;j++)
      { h1=s[j]*s[j];
        ddy[j]=6.0*(y[j+1]-y[j])/h1-2.0*(2.0*dy[j]+dy[j+1])/s[j];
      }
    h1=s[n-2]*s[n-2];
    ddy[n-1]=6.*(y[n-2]-y[n-1])/h1+2.*(2.*dy[n-1]+dy[n-2])/s[n-2];
    g=0.0;
    for (i=0;i<=n-2;i++)
      { h1=0.5*s*(y+y[i+1]);
        h1=h1-s*s*s*(ddy+ddy[i+1])/24.0;
        g=g+h1;
      }
    for (j=0;j<=m-1;j++)
      { h0=t[j];
        while (h0>=x[n-1]) h0=h0-(x[n-1]-x[0]);
        while (h0<x[0]) h0=h0+(x[n-1]-x[0]);
        i=0;
        while (h0>x[i+1]) i=i+1;
        u=h0;
        h1=(x[i+1]-u)/s;
        h0=h1*h1;
        z[j]=(3.0*h0-2.0*h0*h1)*y;
        z[j]=z[j]+s*(h0-h0*h1)*dy;
        dz[j]=6.0*(h0-h1)*y/s;
        dz[j]=dz[j]+(3.0*h0-2.0*h1)*dy;
        ddz[j]=(6.0-12.0*h1)*y/(s*s);
        ddz[j]=ddz[j]+(2.0-6.0*h1)*dy/s;
        h1=(u-x)/s;
        h0=h1*h1;
        z[j]=z[j]+(3.0*h0-2.0*h0*h1)*y[i+1];
        z[j]=z[j]-s*(h0-h0*h1)*dy[i+1];
        dz[j]=dz[j]-6.0*(h0-h1)*y[i+1]/s;
        dz[j]=dz[j]+(3.0*h0-2.0*h1)*dy[i+1];
        ddz[j]=ddz[j]+(6.0-12.0*h1)*y[i+1]/(s*s);
        ddz[j]=ddz[j]-(2.0-6.0*h1)*dy[i+1]/s;
      }
      
    free(s);
    free(dy);
    free(ddy);
    free(dz);
    free(ddz);
    return(g);
  } :'( :'( :'( :'( :'( :'(
您需要登录后才可以回帖 登录 | 注册

本版积分规则

Archiver|小黑屋|联系我们|仿真互动网 ( 京ICP备15048925号-7 )

GMT+8, 2024-5-6 05:10 , Processed in 0.050106 second(s), 13 queries , Gzip On, MemCache On.

Powered by Discuz! X3.5 Licensed

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表