global.h

数值分析的曲线插值算法

#include <iostream>
#include <math.h>
using namespace std;

#define N 13
#define delta 0.5
FILE *file;
float X[13] = {0.00,4.74,9.50,19.00,38.00,57.00,76.00,
95.00,114.00,133.00,152.00,171.00,190.00};
float Y[13] = {0.00,5.32,8.10,11.97,16.15,17.10,16.34,
14.63,12.16,9.69,7.03,3.99,0.00};
float R[(int)(190/delta)] = {0};

float M[13][13] = {0};
float B[13] = {0};

void head_end(float x1,float x2,float *head,float *end)//确定本段取值的起末参数
{
	if(x1-floor(x1)>0.5)*head = ceil(x1);  
	else
		if(x1 == floor(x1))
			*head = x1;
		else 
			*head = floor(x1) + 0.5;
	if(x2-floor(x2)>0.5)*end = ceil(x2)-0.5;
	else 
		if(x2 == floor(x2))
			*end = x2 - 0.5;
		else 
			*end = floor(x2);
}

int linear(float *x,float *y,float *r,int s)  //分段线性插值
{
	float head,end;
	head_end(x[0],x[1],&head,&end);
	float k = (y[1] - y[0])/(x[1] - x[0]);
	float t = head;
	int i=0;
	while(t <= end)	{
		r[i] = y[0]+k*(t - x[0]);
		//		cout<<(s+i)*0.5<<"->"<<r[i]<<"\n";
		t += 0.5;
		i++;
	}
	return s + i;
}

void line(float *x,float *y, float *r)
{
	file = fopen("linear.txt","w");
	int i,j;int s=0;
	for(i=0;i<N-1;i++)
		s = linear(&x[i],&y[i],&r[s],s);
	for(i=0;i<(int)(190/delta);i++)
	{
		if(i%4 == 0)
		{
			fprintf(file,"\n\nx");
			for(j=0;j<4;j++)
				fprintf(file,"\t%f",(i+j)*0.5);
			fprintf(file,"\ny");
		}
		fprintf(file,"\t%f",r[i]);
	}
	fclose(file);
}

int square(float *x,float *y,float *r,int s)
{
	float head,end;
	head_end(x[0],x[2],&head,&end);
	float k1 = (y[1]-y[0])/(x[1]-x[0]);
	float k2 = (y[2]-y[1])/(x[2]-x[1]);
	k2 = (k2 - k1)/(x[2] - x[0]);
	float t = head;
	int i = 0;
	float temp;
	while(t<=end)
	{
		r[i] = y[0] + k1*(t-x[0]) + k2*(t-x[0])*(t-x[1]);
// 		if(s == 0)
// 			r[i] = temp;
// 		else
// 			r[i] = (temp + r[i])/2;
		t += 0.5;
		i++;
	}
	int s_m = s + i;
	/*	head_end(x[1],x[2],&head,&end);
	t = head;
	while(t<=end)
	{
	r[i] = y[0] + k1*(t-x[0]) + k2*(t-x[0])*(t-x[1]);
	t += 0.5;
	i++;
	}
	*/
	return s_m;
}

void sqr(float *x,float *y, float *r)
{
	file = fopen("sqr.txt","w");
	int i,j;int s=0;
	for(i=0;i<N-2;i++)
		s = square(&x[i*2],&y[i*2],&r[s],s);
	for(i=0;i<(int)(190/delta);i++)
	{
		if(i%4 == 0)
		{
			fprintf(file,"\n\nx");
			for(j=0;j<4;j++)
				fprintf(file,"\t%f",(i+j)*0.5);
			fprintf(file,"\ny");
		}
		fprintf(file,"\t%f",r[i]);
	}
	fclose(file);
}

int cubic(float *x,float *y,float *r,int s)
{
	float head,end;
	head_end(x[0],x[3],&head,&end);
	float k1 = (y[1]-y[0])/(x[1]-x[0]);
	float k2 = (y[2]-y[1])/(x[2]-x[1]);
	float k3 = (y[3]-y[2])/(x[3]-x[2]);

	float k4 = (k2 - k1)/(x[2] - x[0]);
	k3 = (k3 - k2)/(x[3] - x[1]);
	k3 = (k3 - k4)/(x[3] - x[0]);

	float t = head;
	int i = 0;
	while(t<=end)
	{
		r[i] = y[0] + k1*(t-x[0]) + k4*(t-x[0])*(t-x[1]) + k3*(t-x[0])*(t-x[1])*(t-x[2]);
		t += 0.5;
		i++;
	}
	return s+i;
}

void cub(float *x,float *y,float *r)
{
	file = fopen("cub.txt","w");
	int i,j;int s=0;
	for(i=0;i<N-3;i++)
		s = cubic(&x[i*3],&y[i*3],&r[s],s);
	for(i=0;i<(int)(190/delta);i++)
	{
		if(i%4 == 0)
		{
			fprintf(file,"\n\nx");
			for(j=0;j<4;j++)
				fprintf(file,"\t%f",(i+j)*0.5);
			fprintf(file,"\ny");
		}
		fprintf(file,"\t%f",r[i]);
	}
	fclose(file);
}

void initM(void)
{
	int i,j;
	M[0][0] = 2;M[N-1][N-1] = 2;
	for(i=1;i<N-1;i++)
	{
		M[i][i] = 2;
		M[i][i+1] = (X[i+1]-X[i])/(X[i+1]-X[i-1]);
		M[i][i-1] = 1-M[i][i+1];
		B[i] = 6/(X[i+1]-X[i-1])*((Y[i+1]-Y[i])/(X[i+1]-X[i])-(Y[i]-Y[i-1])/(X[i]-X[i-1]));
	}
/*
	for(j=0;j<N;j++)
	{
		for(i=0;i<N;i++)
			cout<<M[j][i]<<" ";
		cout<<B[j]<<endl;
	}
//*/
}

void GetM()
{
	int i,j;
	float k;
	for(i=0;i<N-1;i++)
	{
		k = M[i+1][i]/M[i][i];
		M[i+1][i] = 0;
		M[i+1][i+1] -= k*M[i][i+1];
		B[i+1] -= k*B[i];
	}
	/*
	for(j=0;j<N;j++)
	{
	for(i=0;i<N;i++)
	cout<<M[j][i]<<" ";
	cout<<B[j]<<endl;
	}
	cout<<endl;
	//*/
	for(i=N-1;i>1;i--)
	{
		k = M[i-1][i]/M[i][i];
		M[i-1][i] = 0;
		B[i-1] -= k*B[i];
	}
	/*
	for(j=0;j<N;j++)
	{
	for(i=0;i<N;i++)
	cout<<M[j][i]<<" ";
	cout<<B[j]<<endl;
	}
	cout<<endl;
	//*/
	for(i=0;i<N;i++)
	{
		B[i] = B[i]/M[i][i];
		M[i][i] = 1;
	}
	/*
	for(j=0;j<N;j++)
	{
	for(i=0;i<N;i++)
	cout<<M[j][i]<<" ";
	cout<<B[j]<<endl;
	}
	cout<<endl;
	//*/
}

int Spine(float *x,float *y,float *r,int x_i,int s)
{
	float head,end;
	head_end(x[0],x[1],&head,&end);
	float h = x[1]-x[0];
	float k1 = B[x_i]/6/h;
	float k2 = B[x_i+1]/6/h;
	float k3 = y[0]/h-B[x_i]/6*h;
	float k4 = y[1]/h-B[x_i+1]/6*h;
	float t = head;
	int i=0;
	while(t<=end)
	{
		r[i] = k1*(x[1]-t)*(x[1]-t)*(x[1]-t)+
			k2*(t-x[0])*(t-x[0])*(t-x[0])+
			k3*(x[1]-t)+
			k4*(t-x[0]);
		t += 0.5;
		i++;
	}
	return s+i;
}

void Spi(float *x,float *y,float *r)
{
	initM();
	GetM();
	file = fopen("spi.txt","w");
	int i,j;int s=0;
	for(i=0;i<N;i++)
		s = Spine(&x[i],&y[i],&r[s],i,s);
	for(i=0;i<(int)(190/delta);i++)
	{
		if(i%4 == 0)
		{
			fprintf(file,"\n\nx");
			for(j=0;j<4;j++)
				fprintf(file,"\t%f",(i+j)*0.5);
			fprintf(file,"\ny");
		}
		fprintf(file,"\t%f",r[i]);
	}
	fclose(file);
}