interpoints.cpp

vc下实现的分段线性插值、二次多项式插值、三次多项式插值、三次样条插值

#include "InterPoints.h"



//分段线性插值法
BOOL InterLinear(double * pd_X, double * pd_Y, int N, double * pd_resultY)
{
	int n = 1;	//插值阶数
	double * pd_Coe = new double[n+1];	//系数个数为阶数n+1
	int i,k;	//i记录X角标,k记录插值y的位置
	
	i = 0;
	k = 0;
	LinearCoe(pd_X, pd_Y,pd_Coe,i);

	while (k <= pd_X[N-1])
	{	
		if (k > pd_X[i+1])
		{
			i++;
			LinearCoe(pd_X, pd_Y,pd_Coe,i);
		}		
		if (k == pd_X[i+1])
		{
			pd_resultY[k] = pd_Y[i+1];
			k++;
			continue;
		}
		pd_resultY[k] = pd_Coe[0]+pd_Coe[1]*k;
		k++;
	}
	
	return TRUE;
}

//分段二次多项式插值法
BOOL InterQuadrate(double * pd_X, double * pd_Y, int N, double * pd_resultY)
{
	int n = 2;	//插值阶数
	double * pd_Coe = new double[n+1];	//系数个数为阶数n+1
	int i,k;	//i记录X角标,k记录插值y的位置
	double h;	//记录节点步长
	
	i = 1;
	k = 0;
	h = pd_X[2] - pd_X[1];
	QuadrateCoe(pd_X, pd_Y,pd_Coe,i);
	
	while (k <= pd_X[i] + h/2)
	{// x <= x1+h/2的情况	
		pd_resultY[k] = pd_Coe[0]*k*k + pd_Coe[1]*k + pd_Coe[2];
		k++;
	}
	while (k <= pd_X[N-3] + (pd_X[N-2] - pd_X[N-3])/2 )
	{// xj - h/2 < x <= xj + h/2的情况	
		if (k > pd_X[i] + h/2)
		{
			i++;
			QuadrateCoe(pd_X, pd_Y,pd_Coe,i);
			h = pd_X[i+1] - pd_X[i];
		}		
		if (k == pd_X[i])
		{
			pd_resultY[k] = pd_Y[i];
			k++;
			continue;
		}
		pd_resultY[k] = pd_Coe[0]*k*k + pd_Coe[1]*k + pd_Coe[2];
		k++;
	}	
	while (k <= pd_X[N-1])
	{// x > x(n-1) - h/2的情况	
		QuadrateCoe(pd_X, pd_Y,pd_Coe,i);
		if (k == pd_X[i])
		{
			pd_resultY[k] = pd_Y[i];
			k++;
			continue;
		}
		pd_resultY[k] = pd_Coe[0]*k*k + pd_Coe[1]*k + pd_Coe[2];
		if (k == pd_X[N-1])
		{
			pd_resultY[k] = pd_Y[N -1];
		}
		k++;
	}	
	return TRUE;
}

//分段三次多项式插值法
BOOL InterThrice(double * pd_X, double * pd_Y, int N, double * pd_resultY)
{
	int n = 3;	//插值阶数
	double * pd_Coe = new double[n+1];	//系数个数为阶数n+1
	int i,k;	//i记录X角标,k记录插值y的位置
	double h;	//记录节点步长
	
	i = 2;
	k = 0;
	h = pd_X[3] - pd_X[2];
	ThriceCoe(pd_X, pd_Y,pd_Coe,i);
	
	while (k <= pd_X[i] + h/2)
	{// x <= x2+h/2的情况	
		pd_resultY[k] = pd_Coe[0]*k*k*k + pd_Coe[1]*k*k + pd_Coe[2]*k + pd_Coe[3];
		k++;
	}
	while (k <= pd_X[N-4] + (pd_X[N-3] - pd_X[N-4])/2 )
	{// xj - h/2 < x <= xj + h/2的情况	
		if (k > pd_X[i] + h/2)
		{
			i++;
			ThriceCoe(pd_X, pd_Y,pd_Coe,i);
			h = pd_X[i+1] - pd_X[i];
		}		
		if (k == pd_X[i-1])
		{
			pd_resultY[k] = pd_Y[i-1];
			k++;
			continue;
		}
		if (k == pd_X[i])
		{
			pd_resultY[k] = pd_Y[i];
			k++;
			continue;
		}
		pd_resultY[k] = pd_Coe[0]*k*k*k + pd_Coe[1]*k*k + pd_Coe[2]*k + pd_Coe[3];
		k++;
	}	
	int ll = 0;
	while (k <= pd_X[N-1])
	{// x > x(n-1) - h/2的情况	
		ThriceCoe(pd_X, pd_Y,pd_Coe,i);
		if (k == pd_X[i])
		{
			pd_resultY[k] = pd_Y[i];
			k++;
			continue;
		}
		pd_resultY[k] = pd_Coe[0]*k*k*k + pd_Coe[1]*k*k + pd_Coe[2]*k + pd_Coe[3];
		if (k == pd_X[N-1])
		{
			pd_resultY[k] = pd_Y[N -1];
		}
		k++;
	}	
	return TRUE;
}

//三次样条插值法
BOOL Spline(double * pd_X, double * pd_Y, int N, double * pd_resultY)
{
	int n = N - 1;
	double * pd_h = new double[n];	//节点步长数组h
	double * pd_M = new double[N];		//系数数组M
	int i;

	for(i = 0; i < n; i++)	//计算节点步长数组h
	{
		pd_h[i] = pd_X[i + 1] - pd_X[i];
	}	
// 	ShowVector(pd_h,n);

	SplineCoe(pd_X,pd_Y,N,pd_h,pd_M);	
// 	ShowVector(pd_M,N);
 //	ShowVector(pd_h,n);
	
	i = 1;
	int k = 0;

	pd_resultY[k] = pd_Y[0];
	k++;

	while (k <= pd_X[N-1])
	{	
		if (k > pd_X[i])
		{
			i++;
		}		
		if (k == pd_X[i])
		{
			pd_resultY[k] = pd_Y[i];
			k++;
			continue;
		}

		if (k > 171)
		{
			double temp = pow(2,3)/5;
		}

//
		double xi = pd_X[i] - k;
		double xi1 = k - pd_X[i-1];
		double mi = pd_M[i] / (6*pd_h[i-1]);
		double mi1 = pd_M[i-1] / (6*pd_h[i-1]);
		double yi = pd_Y[i]/pd_h[i-1];
		double yi1 = pd_Y[i-1]/pd_h[i-1];
//
//		double tp = xi*xi*xi*mi1 + mi*xi1*xi1*xi1 + (yi1 - mi1) * xi + (yi-mi) * xi1;
//		double tp1 = pow(xi,3)*mi1 + mi*pow(xi1,3) + (yi1 - (mi1*pd_h[i-1]*pd_h[i-1])) * xi + (yi-(mi*pd_h[i-1]*pd_h[i-1])) * xi1;
//
//
//		double x1 = (pow((pd_X[i] - k),3) * pd_M[i-1]) / (6*pd_h[i-1]);
//		double x2 = (pow((k - pd_X[i-1]),3) * pd_M[i]) / (6*pd_h[i-1]);
//		double x3 =((pd_Y[i-1]/pd_h[i-1]) - (pd_M[i-1]*pd_h[i-1]/6)) * (pd_X[i] - k);
//		double x4 =	((pd_Y[i]/pd_h[i-1]) - (pd_M[i]*pd_h[i-1]/6)) * (k - pd_X[i-1]);
// 		double sum = x1+x2+x3-x4;
//		pd_resultY[k] = ((pow((pd_X[i] - k),3) * pd_M[i-1]) / (6*pd_h[i-1])) 
//						+ ((pow((k - pd_X[i-1]),3) * pd_M[i]) / (6*pd_h[i-1])) 
//						+ ((pd_Y[i-1]/pd_h[i-1]) - (pd_M[i-1]*pd_h[i-1]/6)) * (pd_X[i] - k)
//						+ ((pd_Y[i]/pd_h[i-1]) - (pd_M[i]*pd_h[i-1]/6)) * (k - pd_X[i-1]);

//		double c1 = (pd_Y[i] - pd_Y[i-1])/pd_h[i-1] + (pd_M[i-1]-pd_M[i]) * pd_h[i-1]/6;
//		double c2 = pd_Y[i] - ((pd_M[i] * pd_h[i-1] * pd_h[i-1]) /6) - c1*pd_X[i];
//
//		pd_resultY[k] = ((pow((pd_X[i] - k),3) * pd_M[i-1]) / (6*pd_h[i-1])) 
//			+ ((pow((k - pd_X[i-1]),3) * pd_M[i]) / (6*pd_h[i-1]))
//			+ c1*k + c2;

//		pd_resultY[k] = ((pow((pd_X[i] - k),3) * pd_M[i-1]) / (6*pd_h[i-1])) 
//			+ ((pow((k - pd_X[i-1]),3) * pd_M[i]) / (6*pd_h[i-1])) 
//			+ ((pd_Y[i-1]/pd_h[i-1]) - (pd_M[i-1]*pd_h[i-1]/6)) * (pd_X[i] - k)
//			+ ((pd_Y[i]/pd_h[i-1]) - (pd_M[i]*pd_h[i-1]/6)) * (k - pd_X[i-1]);
		

//		pd_resultY[k] = sum;


		pd_resultY[k] = pd_Y[i-1] + (pd_Y[i] - pd_Y[i-1])*(k - pd_X[i-1])/pd_h[i-1]
						- ((pd_M[i] / (6*pd_h[i-1]))/6 + (pd_M[i-1] / (6*pd_h[i-1]))/3) * pd_h[i-1]*(k - pd_X[i-1]) 
						+ pow((k - pd_X[i-1]),2) * (pd_M[i-1] / (6*pd_h[i-1]))/2
						+ ((pd_M[i] / (6*pd_h[i-1])) - (pd_M[i-1] / (6*pd_h[i-1]))) * pow((k - pd_X[i-1]),3) / (6 * pd_h[i-1]);

	

//		pd_resultY[k] = pd_Y[i-1] + (pd_Y[i] - pd_Y[i-1])*(k - pd_X[i-1])/pd_h[i-1]
//			- ((pd_M[i] / (6*pd_h[i-1]))/6 + (pd_M[i-1] / (6*pd_h[i-1]))/3) * pd_h[i-1]*(k - pd_X[i-1]) 
//			+ pow((k - pd_X[i-1]),2) * (pd_M[i-1] / (6*pd_h[i-1]))/2
//			+ ((pd_M[i] / (6*pd_h[i-1])) - (pd_M[i-1] / (6*pd_h[i-1]))) * pow((k - pd_X[i-1]),3) / (6 * pd_h[i-1]);

		if (k == 189)
		{
			double temp = pow(2,3)/5;
		}
		
		k++;
	}
	

	return TRUE;

}

void main()
{
	int N = 13;	//原始数据个数
	int Num = 191;	//插值后Y的数据，X从0到190，增量为1
	double Arrayd_X[13] = {0,4.74,9.5,19,38,57,76,95,114,133,152,171,190};	//原始数据X
	double Arrayd_Y[13] = {0,5.32,8.1,11.97,16.15,17.1,16.34,14.63,12.16,9.69,7.03,3.99,0};	//原始数据Y

//	int N = 4;	//原始数据个数
//	int Num = 20;	//插值后Y的数据，X从0到190，增量为1
//	double Arrayd_X[4] = {0,4.5,9,13.5};	//原始数据X
//	double Arrayd_Y[4] = {0,5,8,13};	//原始数据Y
	

	double * pd_resultY = new double[Num];
	
	Spline(Arrayd_X,Arrayd_Y,N,pd_resultY);
//	InterThrice(Arrayd_X,Arrayd_Y,N,pd_resultY);
//	InterQuadrate(Arrayd_X,Arrayd_Y,N,pd_resultY);
//	InterLinear(Arrayd_X,Arrayd_Y,N,pd_resultY);
	for(int i = 0; i < Num; i++)
	{
		cout<<pd_resultY[i]<<endl;
	}

//	cout<<pow(2,3);
}