三次样条插值的C语言实现.rar

#include<stdio.h> #include<math.h> int main(void){ double x[19]={0.52,3.1,8.0,17.95,28.65,39.62,50.65,78,104.6,156.6,208.6,260.7,312.5,364.4,416.3,468,494,507,520}; double y[19]={5.28794,9.4,13.84,20.2,24.9,28.44,31.1,35,36.5,36.6,34.6,31.0,26.34,20.9,14.8,7.8,3.7,1.5,0.2}; double t[12]={2,4,6,12,16,30,60,110,180,280,400,515}; double h[19],u[19],l[19],g[19]; double M[19],n[19]; void tri(double a[],double d[],double c[],double b[],double M[]); int i,j; double result; for(i=1;i<=18;i++) h[i]=x[i]-x[i-1]; for(i=1;i<=17;i++) u[i]=h[i]/(h[i]+h[i+1]); for(i=1;i<=17;i++) l[i]=1-u[i]; for(i=1;i<=17;i++) g[i]=6.0/(h[i]+h[i+1])*((y[i+1]-y[i])/h[i+1]-(y[i]-y[i-1])/h[i]); g[0]=6.0/h[1]*((y[1]-y[0])/h[1]-1.86548); g[18]=6.0/h[18]*((-0.046115)-(y[18]-y[17])/h[18]); u[18]=1; l[0]=1; for(i=0;i<=18;i++) n[i]=2; tri(u,n,l,g,M); for(i=0;i<12;i++) { for(j=0;j<18;j++) if(t[i]>=x[i]&&t[i]<=x[i+1])break; result=M[i]*pow((x[i+1]-t[i]),3)/(6*h[i+1])+M[i+1]*pow((t[i]-x[i]),3)/(6*h[i+1])+(y[i]-M[i]*h[i+1]*h[i+1]/6)*(x[i+1]-t[i])/h[i+1]+(y[i+1]-M[i+1]*h[i+1]*h[i+1]/6)*(t[i]-x[i])/h[i+1]; printf("%f ",result); } } void tri(double a[],double d[],double c[],double b[],double M[]) /*使用追赶法求解方程组Mi*/ { int i,n; double xmult; for(i=1;i<=18;i++) { xmult=a[i]/d[i-1]; d[i]-=xmult*c[i-1]; b[i]-=xmult*b[i-1]; } M[18]=b[18]/d[18]; for(i=17;i>=0;i--) { M[i]=(b[i]-c[i]*M[i+1])/d[i]; } }