3107002005_5th_spline.cpp

曲线的三次样条插值

// 3107002005_5th_spline.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include <iostream.h>
#include<fstream.h>
#include<math.h>

#define n 5;
double x[21]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
double y[21]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
double a[21]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
double b[21]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
double c[21]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
double d[21]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
double u[21]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
double l[21]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
double yy[21]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
double m[21]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
double y_0=50/26/26;
double d0=5.25;
double dn=d0;
double un=1;
double v0=1;


int spline(int q=0)
{
	double h[21];
	for(int i=0;i<=5;i++)
	{
		m[i]=y[i];
	}
	for(int k=1;k<=2;k++)
	{
		for(i=5;i>=k;i--)
		{
			m[i]=(m[i]-m[i-1])/(x[i]-x[i-1]);
		}
	}
	h[1]=x[1]-x[0];
	for(i=1;i<=5-1;i++)
	{
		h[i+1]=x[i+1]-x[i];
		c[i]=h[i+1]/(h[i]+h[i+1]);
		a[i]=1-c[i];
		b[i]=2;
		m[i]=6*m[i+1];
	}
	m[0]=d0;
	m[5]=dn;
	c[0]=v0;
	b[0]=2;
	a[5]=un;
	b[5]=2;
	for(i=0;i<=5;i++)
		d[i]=m[i];
	return 0;
}
double compute(double q)
{
	int i=q/20;
	double h=0.4;
	double X1=x[i+1]-(-1+q*2/100);
	double X2=-1+q*2/100-x[i];
	double Y1;double Y2;double Y3;
	Y1=m[i];Y2=m[i+1];
	Y3=Y1*X1*X1*X1*1/6+Y2*X2*X2/6+(y[i]-Y1*h*h/6)*X1+(y[i+1]-Y2*h*h/6)*X2;
	if(q<50)
		a[((int)q)]=Y3/h;
	else
		a[(int)q-79]=Y3/h;
	return Y3/h;
}




void tss(int nn)
{
	u[0]=b[0];
	yy[0]=d[0];
	for(int k=1;k<=nn;k++)
	{
		l[k]=a[k]/u[k-1];
		u[k]=b[k]-l[k]*c[k-1];
		yy[k]=d[k]-l[k]*yy[k-1];
	}
	m[nn]=yy[nn]/u[nn];
	for(k=nn-1;k>=0;k--)
	{
		m[k]=(yy[k]-c[k]*m[k+1])/u[k];

	}
}



int main(int argc, char* argv[])
{
	cout<<((0.04+0.125)/0.4-0.78)*15;
	ofstream outfile1;
	double error=0;
	
	double y_n=-50/26/26;


	outfile1.open("spline插值输出.txt");
	cout<<endl<<"函数在各个插值点的函数值为："<<endl;
	double i;
	for(int i1=1;i1<=10;i1++)
	{
		i=i1;
		y[i1]=1/(1+25*(-1+i*2/100)*(-1+i*2/100));
		x[i1]=-1+i*2/5;
		cout<<y[i1]<<' ';
	}
	for( i1=0;i1<=5;i1++)
	{
		i=i1;
		y[i1]=1/(1+25*(-1+i*2/5)*(-1+i*2/5 ));
		x[i1]=-1+i*2/5;
	//	cout<<y[i1]<<' ';
	}
	spline();
	tss(5);
	cout<<endl<<"样条插值得到的系数矩阵为"<<endl<<"{";
	for(i1=0;i1<=5;i1++)
		cout<<m[i1]<<' ';
	cout<<"}"<<endl;
	cout<<endl<<"插值函数得到的各个数据点的函数值为："<<endl;
	for(double ii=1;ii<=10;ii++)
	{
		i1=ii;//cout<<compute(1);
		compute(ii);
		cout<<a[i1]<<' ';
	}
	cout<<endl;
	for(i=99;i>=90;i--)
	{
		i1=i;
		compute(i);
		cout<<a[i1-79]<<' ';
	}
	cout<<endl;
	//cout<<compute(90);



	return 0;
}