s.cpp

计算三次样条插值函数

#include<iostream.h>

 
void main()
{
	int i;
	long double b[21],x[21],f[21],l,m[21],s[81],d[21],a[20],c[20];
	for(i=0;i<21;i++){
		x[i]=-4.0+i*0.4;
		f[i]=1.0/(x[i]*x[i]*x[i]*x[i]+1.0);
	}
	for(i=0;i<20;i++){
		a[i]=1.0;
		c[i]=1.0;
	}
	d[0]=2.0;
	d[20]=2.0;
	for(i=1;i<20;i++){
		d[i]=4.0;
	}
	b[0]=15.0*((f[1]-f[0])*5.0/2.0-256.0/(257.0*257.0));
	b[20]=15.0*(-256.0/(257.0*257.0)-(f[20]-f[19])*5.0/2.0);
	for(i=1;i<20;i++){
		b[i]=75.0/2.0*(f[i+1]-2.0*f[i]+f[i-1]);
	}
	for(i=0;i<20;i++){
		l=a[i]/d[i];
		d[i+1]=d[i+1]-l*c[i];
		b[i+1]=b[i+1]-l*b[i];
	}
	m[20]=b[20]/d[20];
	for(i=19;i>=0;i--){
		m[i]=(b[i]-c[i]*m[i+1])/d[i];
	}
	for(i=0;i<20;i++){
		s[4*i]=m[i]*5.0/12.0*(x[i+1]-x[i])*(x[i+1]-x[i])*(x[i+1]-x[i])+(f[i]*5.0/2.0-m[i]/15.0)*(x[i+1]-x[i]);
		s[4*i+1]=m[i]*5.0/12.0*(x[i+1]-x[i]-0.1)*(x[i+1]-x[i]-0.1)*(x[i+1]-x[i]-0.1)+m[i+1]*5.0/12.0*(x[i]+0.1-x[i])*(x[i+1]-x[i]+0.1)*(x[i+1]-x[i]+0.1)+(f[i+1]*5.0/2.0-m[i+1]/15.0)*0.1+(f[i]*5.0/2.0-m[i]/15.0)*(x[i+1]-x[i]-0.1);
		s[4*i+2]=m[i]*5.0/12.0*(x[i+1]-x[i]-0.2)*(x[i+1]-x[i]-0.2)*(x[i+1]-x[i]-0.2)+m[i+1]*5.0/12.0*(x[i]+0.2-x[i])*(x[i+1]-x[i]+0.2)*(x[i+1]-x[i]+0.2)+(f[i+1]*5.0/2.0-m[i+1]/15.0)*0.2+(f[i]*5.0/2.0-m[i]/15.0)*(x[i+1]-x[i]-0.2);
		s[4*i+3]=m[i]*5.0/12.0*(x[i+1]-x[i]-0.3)*(x[i+1]-x[i]-0.3)*(x[i+1]-x[i]-0.3)+m[i+1]*5.0/12.0*(x[i]+0.3-x[i])*(x[i+1]-x[i]+0.3)*(x[i+1]-x[i]+0.3)+(f[i+1]*5.0/2.0-m[i+1]/15.0)*0.3+(f[i]*5.0/2.0-m[i]/15.0)*(x[i+1]-x[i]-0.3);
	}
	s[80]=m[20]*5.0/12.0*(x[20]-x[19])*(x[20]-x[19])*(x[20]-x[19])+(f[20]*5.0/2.0-m[20]/15.0)*(x[20]-x[19]);
	l=-4.0;
	for(i=0;i<81;i++){
		cout<<l<<'\t'<<s[i]<<endl;
		l=l+0.1;
	}
}