interpoints.h

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

#include "math.h"
#include "stdio.h"
#include <iomanip.h>
#include "iostream.h"
#include "GaussRemove.h"

#define BOOL int
#define FALSE 0
#define TRUE 1

//计算分段线性插值的系数
void LinearCoe(double * pd_X, double * pd_Y,double * pd_Coe, int i)
{
	pd_Coe[0] = (-pd_X[i]*pd_Y[i+1]+pd_X[i+1]*pd_Y[i])/(pd_X[i+1]-pd_X[i]);
	pd_Coe[1] = (pd_Y[i+1]-pd_Y[i])/(pd_X[i+1]-pd_X[i]);
}

//计算分段二次多项式插值的系数
void QuadrateCoe(double * pd_X, double * pd_Y,double * pd_Coe, int i)
{
	int n = 2;
	double * pd_MatrixA = new double[(n+1)*(n+1)];
	double * pd_VectorB = new double[n+1];

	pd_VectorB[0] = pd_Y[i-1];
	pd_VectorB[1] = pd_Y[i];
	pd_VectorB[2] = pd_Y[i+1];

	for(int j = 0; j < n+1; j++)
	{
		pd_MatrixA[j*3] = pd_X[i+j-1]*pd_X[i+j-1];
		pd_MatrixA[j*3 + 1] = pd_X[i+j-1];
		pd_MatrixA[j*3 + 2] = 1;
	}
	GaussRemove(pd_MatrixA,pd_VectorB,pd_Coe,n+1);
}

//计算分段三次多项式插值的系数
void ThriceCoe(double * pd_X, double * pd_Y,double * pd_Coe, int i)
{
	int n = 3;
	double * pd_MatrixA = new double[(n+1)*(n+1)];
	double * pd_VectorB = new double[n+1];
	
	pd_VectorB[0] = pd_Y[i-2];
	pd_VectorB[1] = pd_Y[i-1];
	pd_VectorB[2] = pd_Y[i];
	pd_VectorB[3] = pd_Y[i+1];
	
	for(int j = 0; j < n+1; j++)
	{
		pd_MatrixA[j*4] = pd_X[i+j-2]*pd_X[i+j-2]*pd_X[i+j-2];
		pd_MatrixA[j*4 + 1] = pd_X[i+j-2]*pd_X[i+j-2];
		pd_MatrixA[j*4 + 2] = pd_X[i+j-2];
		pd_MatrixA[j*4 + 3] = 1;		
	}
	GaussRemove(pd_MatrixA,pd_VectorB,pd_Coe,n+1);
}

//计算三次样条插值函数的系数组M
void SplineCoe(double * pd_X, double * pd_Y, int N, double * pd_h, double * pd_M)
{
	int n = N -1;
	double * pd_Alpha = new double[n];	//系数数组Alpha
	double * pd_Beta = new double[n];	//系数数组Beta
	double * pd_Gama = new double[n];	//系数数组Gama
	int i,j;

	for(i = 0; i < n; i++)		//计算Alpha,Beta,Gama
	{
		pd_Alpha[i] = pd_h[0]/(pd_h[0] + pd_h[i]);
		pd_Gama[i] = 1 - pd_Alpha[i];
		pd_Beta[i] = 6 * (((pd_Y[1]-pd_Y[0])/pd_h[0]) - ((pd_Y[i+1]-pd_Y[i])/pd_h[i])) / (pd_h[0] + pd_h[i]);
	}
	
	double * pd_MatrixA = new double[n*n];
	
	for(i = 0; i < n; i++)
	{
		for(j = 0; j < n; j++)
		{
			if (i == j)
			{
				pd_MatrixA[i*n+j] = 2;
			}
			else if (j == (i+1))
			{
				pd_MatrixA[i*n+j] = pd_Alpha[i];			
			}
			else if (j == (i-1))
			{
				pd_MatrixA[i*n+j] = pd_Gama[i];			
			}
			else
				pd_MatrixA[i*n+j] = 0;
		}
	}
	
	pd_MatrixA[n-1] = pd_Gama[0];
	pd_MatrixA[(n-1)*n] = pd_Alpha[n-1];	

//	ShowMatrix(pd_MatrixA,n,n);
	GaussRemove(pd_MatrixA,pd_Beta,&(pd_M[1]),N-1);
	pd_M[0] = pd_M[N-1];

// 	ShowVector(pd_Beta,n);

	delete pd_Alpha;
	delete pd_Beta;
	delete pd_Gama;
	delete pd_MatrixA;

	
	
}