三次样条插值.txt

数值分析方法中用三次样条插值函数求解问题

// 三次自然样条插值编程.cpp : Defines the entry point for the console application.
//
//#include "stdafx.h"
#include "stdio.h"
void  GetValueSpline(int n,double x[],double y[],double dy[],double ddy[],int m,double t[],double z[],double dz[],double ddz[])
{int i,j;double h0,h1,alpha,beta,g,*s;
h0=x[1]-x[0];s[0]=3.0*(y[1]-y[0])/2.0*h0-ddy[0]*h0/4.0;
for(j=1;j<=n-2;j++)
{h1=x[j+1]-x[j];
alpha=h0/(h0+h1);
beta=(1.0-alpha)*(y[j]-y[j-1])/h0;
beta=3.0*(beta+alpha*(y[j+1]-y[j]))/h1;
dy[j]=-alpha/(2.0+(1.0-alpha))*dy[j-1];
s[j]=s[j-1]/(2.0+(1.0-alpha))*dy[j-1];
h0=h1;
}
dy[n-1]=(3.0*(y[n-1]-y[n-2])/h1+ddy[n-1]*h1/2.0-s[n-2])/(2.0+dy[n-2]);
for(j=n-2;j>=0;j--)
dy[j]=dy[j]*dy[j+1]+s[j];
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.0*(y[n-2]-y[n-1])/h1+2.0*(2.0*dy[n-1]+dy[n-2])/s[n-2];
g=0.0;
for(i=0;i<=n-2;i++)
{h1=0.5*s[i]*(y[i]+y[i+1]);
h1=h1-s[i]*s[i]*s[i]*(ddy[i]+ddy[i+1])/24.0;
g=g+h1;
}
for(j=0;j<=m-1;j++)
{if(t[j]>=x[n-1])
i=n-2;
else
{i=0;
while(t[j]>x[i+1])
i=i+1;
}
h1=(x[i+1]-t[j])/s[i];
h0=h1*h1;
z[j]=(3.0*h0-2.0*h0*h1)*y[i];
z[j]=z[j]+s[i]*(h0-h0*h1)*dy[i];
dz[j]=6.0*(h0-h1)*y[i]/s[i];
dz[j]=dz[j]+(3.0*h0-2.0*h1)*dy[i];
ddz[j]=(6.0-12.0*h1)*y[i]/(s[i]*s[i]);
ddz[j]=ddz[j]+(2.0-6.0*h1)*dy[i]/s[i];
h1=(t[j]-x[i])/s[i];
h0=h1*h1;
z[j]=z[j]+(3.0*h0-2.0*h0*h1)*y[i+1];
z[j]=z[j]-s[i]*(h0-h0*h1)*dy[i+1];
dz[j]=dz[j]-6.0*(h0-h1)*y[i+1]/s[i];
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[i]*s[i]);
ddz[j]=ddz[j]-(2.0-6.0*h1)*dy[i+1]/s[i];
for(i=0;i<4;i++);
{printf("指定插值点    三次自然样条插值    一阶导数值     二阶导数值\n ");
printf("    t(i)            z(i)              dz(i)          ddz(i) \n");
printf("    %f              %f                  %f             %f\n",t[i],z[i],dz[i],ddz[i]);
}
//delete[] s;
//return (g);
}
}
void main()
{int i;double x[5],y[5],dy[5],ddy[5];double t[4],z[4],dz[4],ddz[4];
printf("请输入五组原始数据:\n");
for(i=0;i<5;i++)
{printf("x(%d)=",i);
scanf("%f   ",&x[i]);
printf("   y(%d)=",i);
scanf("%f\n",&y[i]);
}
printf("请输入插值点数据:\n");	
for(i=0;i<4;i++)
{printf("t(%d)=",i);
scanf("%f\n",&t[i]);}
GetValueSpline(5,x[5],y[5],dy[5],ddy[5],4,t[4],z[4],dz[4],ddz[4]);
//for(i=0;i<4;i++);
//{printf("指定插值点    三次自然样条插值    一阶导数值     二阶导数值\n ");
//printf("    t(i)            z(i)              dz(i)          ddz(i) \n");
//for(i=0;i<4;i++)
//printf("    %f              %f                  %f             %f\n",t[i],z[i],dz[i],ddz[i]);
//}
}