c2三次样条插值.txt

C2三次样条插值． ＭＡＴＬＡＢ源程序

#include<stdio.h>
#include<math.h>
#include<alloc.h>
 
zhuiganfa(n,a,b,c,f,x)
    int n;
    double *a,*b,*c,*f,*x;
{
    int i;
    double *B,*y;
    B=malloc(n*sizeof(double));
    y=malloc(n*sizeof(double));
    B[0]=c[0]/b[0];
    for(i=1;i<n-1;i++)
       B[i]=c[i]/(b[i]-a[i]*B[i-1]);
    y[0]=f[0]/b[0];
    for(i=1;i<n;i++)
       y[i]=(f[i]-a[i]*y[i-1])/(b[i]-a[i]*B[i-1]);
    x[n-1]=y[n-1];
    for(i=n-2;i>=0;i--)
       x[i]=y[i]-B[i]*x[i+1];
    free(B);
    free(y);
}
 
double SplineIntp2(int n,double *X,double *Y,double M0,double Mn,double x)
{
    int i;
    double *a,*b,*c,*g,*M;
    double h0,h1,f0,f1,d0,d1,y;
    a=malloc((n+1)*sizeof(double));
    b=malloc((n+1)*sizeof(double));
    c=malloc((n+1)*sizeof(double));
    g=malloc((n+1)*sizeof(double));
    M=malloc((n+1)*sizeof(double));
    b[0]=b[n]=2;
    a[n]=c[0]=1;
    for(i=1;i<n;i++)
    {
       h0=X[i]-X[i-1]; h1=X[i+1]-X[i];
       f0=(Y[i]-Y[i-1])/h0; f1=(Y[i+1]-Y[i])/h1;
       a[i]=h0/(h0+h1); b[i]=2; c[i]=h1/(h0+h1);
       g[i]=6*(f1-f0)/(h0+h1);
    }
    g[1]=g[1]-a[1]*M0;
    g[n-1]=g[n-1]-c[n-1]*Mn;
    zhuiganfa(n-1,a+1,b+1,c+1,g+1,M+1);
    M[0]=M0;
    M[n]=Mn;
    i=1;
    while(x>X[i]) i++;
    h0=X[i]-X[i-1]; d0=x-X[i-1]; d1=X[i]-x;
    y=M[i-1]*d1*d1*d1/(6*h0)+M[i]*d0*d0*d0/(6*h0)+(Y[i-1]-M[i-1]*h0*h0/6)*d1/h0+(Y[i]-M[i]*        h0*h0/6)*d0/h0;
    free(a); free(b); free(c); free(g); free(M);
    return y;
}
 
main()
{
   double X[4]={0,1,2,3},Y[4]={0,0,0,0},M0=1,Mn=0,x1=0.5,x2=1.5,x3=2.5,y1,y2,y3;
   y1=SplineIntp2(3,X,Y,M0,Mn,x1);
   y2=SplineIntp2(3,X,Y,M0,Mn,x2);
   y3=SplineIntp2(3,X,Y,M0,Mn,x3);
   printf("y(0.5)=%f\n",y1);
   printf("y(1.5)=%f\n",y2);
   printf("y(2.5)=%f\n",y3);
}