yangtiao.cpp

数值分析里的样条插值算法

// yangtiao.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"

#define N 16

int main(int argc, char* argv[])
{
	int i;	
	double x1[N],y1[N],a[N-1],b[N],c[N-1],h[N-1],p[N-1],q[N-1],y1=0.9994,y2=0.4695;	
	double x[N]={0.03491,0.10472,0.17453,0.24435,0.31416,0.38397,0.45379,0.52360,
				0.59341,0.66323,0.73304,0.80285,0.87266,0.94248,1.01229,1.08210},
			y[N]={0.0349,0.1045,0.1736,0.2419,0.3090,0.3746,0.4384,0.5000,
				0.5592,0.6157,0.6691,0.7193,0.7660,0.8090,0.8480,0.8826};

	a[0]=c[N-2]=1;
	for(i=0;i<N-1;i++)
	{
		h[i]=x[i+1]-x[i];
	}
	b[0]=1.0*6/h[0]*((y[1]-y[0])/h[0]-y1);
	b[N-1]=1.0*6/h[N-2]*(y2-(y[N-1]-y[N-2])/h[N-2]);
	for(i=1;i<N-1;i++)
	{
		a[i]=h[i]/(h[i-1]+h[i]);
		b[i]=1.0*6/(h[i-1]+h[i])*((y[i+1]-y[i])/h[i]-(y[i]-y[i-1])/h[i-1]);
		c[i-1]=1-a[i];
	}
	
	p[0]=2;
	y1[0]=b[0]/p[0];
	for(i=0;i<N-1;i++)
	{
		q[i]=a[i]/p[i];
		p[i+1]=2-c[i]*q[i];
		y1[i+1]=(b[i+1]-c[i]*y1[i])/p[i+1];
	}
	x1[N-1]=y1[N-1];
	for(i=N-2;i>=0;i--)
	{
		x1[i]=y1[i]-q[i]*x1[i+1];
	}
	for(i=0;i<N-1;i++)
		printf("x1[%d]=%5.4f\t",i+1,x1[i]);

 

	return 0;
}