三样条插值(第一种边界条件下)

#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 SplineIntp1(int n,double *X,double *Y,double ma,double mb,double x)

    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;

    {

       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);

    }

    h0=X[1]-X[0]; g[0]=6*((Y[1]-Y[0])/h0-ma)/h0;

    h1=X[n]-X[n-1]; g[n]=6*(

mb-

(Y[n]-Y[n-1])/h1)/h1;

    zhuiganfa(n+1,a,b,c,g,M);

    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},ma=1,mb=0,x1=0.5,x2=1.5,x3=2.5,y1,y2,y3;

   y1=SplineIntp1(3,X,Y,ma,mb,x1);

   y2=SplineIntp1(3,X,Y,ma,mb,x2);

   y3=SplineIntp1(3,X,Y,ma,mb,x3);

   printf("y(0.5)=%f\n",y1);

   printf("y(1.5)=%f\n",y2);

   printf("y(2.5)=%f\n",y3);