gaussremove.h

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

#include "math.h"
#include "iostream.h"

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


//显示矩阵
void ShowMatrix(double * P_Matrix,int M, int N)
{
	cout<<endl;
	for(int i = 0; i < M; i++)
	{
		for(int j = 0; j < N; j++)
		{
			
			cout<<P_Matrix[i*N+j]<<"   ";
			if (j == N-1)
			{
				cout<<endl;
			}
		}
	}	
}

//显示向量
void ShowVector(double * P_Vector,int M)
{
	cout<<endl;
	
	for(int i = 0; i < M; i++)
	{
		cout<<P_Vector[i]<<endl;
	}
	
	
}


//对矩阵进行顺序高斯消去
BOOL GaussRe(double * P_Matrix,double * P_MatrixG,double * P_Vector, double * P_VectorG,int M, int N)
{
	double m,temp;
	
	int i,j,k,l;

	int flag = 0;
	
//	if (P_Matrix[0] == 0)
//	{
//		return FALSE;
//	}

	//复制原始矩阵和向量

	for(i = 0; i < M; i++)
	{
		P_VectorG[i] = P_Vector[i];
		for(j = 0; j < N; j++)
		{
			P_MatrixG[i*N+j] = P_Matrix[i*N+j];
		}
	}


	//开始消去法
	for(k = 0; k < M-1; k++)
	{
		temp = P_MatrixG[k];
		flag = k;


		//找出主元素
		for(l = k; l < M; l++ )
		{
			if (fabs(P_MatrixG[l*N+k]) > fabs(temp)) 
			{
				temp = P_MatrixG[l*N+k];
				flag = l;
			}
		}

		//交换行
		if (flag != k)
		{
			for(l = k; l < N; l++)
			{
				//交换系数矩阵元素
				temp = P_MatrixG[flag*N+l];
				P_MatrixG[flag*N+l] = P_MatrixG[k*N+l];
				P_MatrixG[k*N+l] = temp;
			}

			
			//交换向量元素
			temp = P_VectorG[flag];
			P_VectorG[flag] = P_VectorG[k];
			P_VectorG[k] = temp;
		}
		
		
		
		for(i = k+1; i < M; i++)
		{			
			m = P_MatrixG[i*N+k]/P_MatrixG[k*N+k];

			P_MatrixG[i*N+k] = 0;	//消去的元素为0

			
			for(j = k+1; j < N; j++)
			{
				P_MatrixG[i*N+j] = P_MatrixG[i*N+j]-m*P_MatrixG[k*N+j];
			}

			P_VectorG[i]= P_VectorG[i] - m*P_VectorG[k];
		}

	}

	return TRUE;

}


//回代求解
BOOL Calculate(double * P_MatrixG,double * P_VectorG, double * P_VectorX,int M, int N)
{
	int j,k;
	double temp;

	P_VectorX[N-1] = P_VectorG[M-1]/P_MatrixG[(M-1)*N+N-1];

	for(k = N-2; k >= 0; k--)
	{
		temp = 0;
		for(j = k+1; j < N; j++)
			{
				temp += P_MatrixG[k*N+j]*P_VectorX[j];
			}
		P_VectorX[k] = (P_VectorG[k] - temp)/P_MatrixG[k*N+k];
	}
	
	return TRUE;
	
}


BOOL GaussRemove(double * pd_MatrixA,double* pd_VectorB,double * pd_Coe, int n)
{
//	int M,N;
//	int i,j;
	
//	cout<<"Input the row and column of the Matrix:"<<endl;
//
//	cin>>M>>N;

//	N = 3;
//
//	double  = new double[N*N];
//	double  = new double[N];
//	double * pf_VectorX = new double[N];

	
//	for(i = 0; i < M; i++)
//	{
//		cout<<"Input the row "<<i+1<<" number of the Matrix:"<<endl;
//		for(j = 0; j < N; j++)
//		{
//			cin>>pf_MatrixA[i*N+j];
//		}
//	}

	
//	pf_MatrixA[0] = 1;
//	pf_MatrixA[1] = 3.5;
//	pf_MatrixA[2] = 1;
//	pf_MatrixA[3] = 1;
//	pf_MatrixA[4] = 1;
//	pf_MatrixA[5] = 3;
//	pf_MatrixA[6] = 2;
//	pf_MatrixA[7] = 3;
//	pf_MatrixA[8] = 3;
//
//	pf_VectorB[0] = 0;
//	pf_VectorB[1] = 0;
//	pf_VectorB[2] = 0;
//	pf_VectorB[0] = 0.6781;
//	pf_VectorB[1] = 12.1;
//	pf_VectorB[2] = 981;
	
//	for(i = 0; i < N; i++)
//	{
//		for(j = 0; j < N; j++)
//		{
//			if (i == j)
//			{
//				pd_MatrixA[i*N+j] = pd_MatrixA[i*N+j] - 6.3236;
//			}
//		}
//	}



//		cout<<"Input the vector B:";
//		for(i = 0; i < M; i++)
//		{
//			cin>>pf_VectorB[i];
//		}


	
// 	ShowMatrix(pd_MatrixA, n,n);

	double * pd_MatrixA_Gauss = new double[n*n];
	double * pd_VectorB_Gauss = new double[n];

	GaussRe(pd_MatrixA,pd_MatrixA_Gauss,pd_VectorB,pd_VectorB_Gauss,n,n);

	Calculate(pd_MatrixA_Gauss, pd_VectorB_Gauss, pd_Coe, n, n);

	delete pd_MatrixA_Gauss;
	delete pd_VectorB_Gauss;


//	ShowMatrix(pd_MatrixA_Gauss, N,N);

//	ShowVector(pf_VectorB_Gauss,M);
	
//	ShowVector(pd_Coe,N);

	return TRUE;
	
}