main.cpp

计算方法中最小二乘法程序设计的源代码,具有很强的通用性

#include "stdio.h"
#include "math.h"

//#define M 11	//采样个数
//#define N 5		//未知数个数

double	dt[]={	1, 25, 8.4634, 6.3548, 1.5561, 3.6599,
				1, 35, 5.8744, 4.6514, 1.2122, 2.6261,
				1, 45, 4.6508, 3.7096, 1.0355, 2.0555, 
				1, 55, 3.4163, 2.8078, 0.9106, 1.6241, 
				1, 65, 2.8506, 2.2576, 0.8349, 1.3675, 
				1, 75, 2.2308, 1.9725, 0.7325, 1.1514, 
				1, 85, 1.7921, 1.6598, 0.6401, 0.9599, 
				1, 95, 1.5319, 1.4081, 0.6160, 0.8654, 
				1, 105,1.3758, 1.1834, 0.5707, 0.7766, 
				1, 115,1.1659, 1.1148, 0.5270, 0.7277,
				1, 125,0.9435, 0.9263, 0.4800, 0.6331};
//double	B[N];
//double	YJ[M];
//double	A[N][N+1];

void LeastSquare(double * Data, short M, short N)
{
	int	i, j, k;

	//定义变量
	double (*A)[6];
	A = new double[N][6];

	double (*X)[5];
	X = new double[M][5];

	double * B = new double[N];
	double YJ[11];
	
	//取数据
	for(i=0; i<M; i++)
	{
		for(j=0; j<N;j++)
		{
			X[i][j]=Data[i*(N+1)+j];
		}
	}
	//计算增广矩阵元素
	
				
	for(k=0; k<N; k++)
	{
		for(j=0; j<N; j++)
		{
			A[k][j] = 0;
						
			for(i=0; i<M; i++)
				A[k][j] = A[k][j]+X[i][k]*X[i][j];
		}
					
		A[k][N] = 0;
		
		for(i=0; i<M; i++)
			A[k][N] = A[k][N]+Data[i*(N+1)+N]*X[i][k];
	}
				
	//列主元高斯－约当消去法解方程组
	double	c, ftemp;
	int		p;
				
	for(k=0;k<N;k++)
	{
		c=A[k][k];
		p=k;
					
		for(i=k+1; i<N; i++)
		{
			if(fabs(A[i][k]) > c)
			{
				c=A[i][k];
				p=i;
			}
		}
		
		if (fabs(c) < 1.0e-9) return;
					
		for(j=k; j<N+1; j++){
			ftemp=A[k][j];
			A[k][j]=A[p][j];
			A[p][j]=ftemp;
		}
					
		for(i=0; i<N; i++){
			if (i==k) continue;
						
			for(j=k+1; j<N+1; j++)
				A[i][j]=A[i][j]-A[i][k]*A[k][j]/A[k][k];
		}
	}
				
	for(i=0; i<N; i++)
		B[i]=A[i][N]/A[i][i];
				
				
	for(i=0; i<M; i++)
	{
		YJ[i]=0;
		for(j=0; j<N; j++)
			YJ[i]+=B[j]*X[i][j];
	}
				
	p=0;

	return;
}

void main()
{	

	LeastSquare( dt, 11, 5 );

	return;
}