lequations.cpp

一个最小二乘法的拟合数据的小程序

// LEquations.cpp: implementation of the CLEquations class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "LEquations.h"
#include "math.h"

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CLEquations::CLEquations()
{

}

CLEquations::~CLEquations()
{

}



//指定系数和常数的构造函数
CLEquations::CLEquations(const CMatrix& mtxCoef,const CMatrix& mtxConst)
{
	ASSERT(Init(mtxCoef,mtxConst));
}

//初始化
BOOL CLEquations::Init(const CMatrix& mtxCoef,const CMatrix& mtxConst)
{
	if (mtxCoef.GetNumRows()!=mtxConst.GetNumRows()) return FALSE;
	else
	{
		m_mtxCoef=mtxCoef;
		m_mtxConst=mtxConst;
		
		return TRUE;
	}


}

/*

//
//获取属性
//

CMatrix GetCoefMatrix() const;//获取系数矩阵
CMatrix GetConstMatrix() const;//获取常数矩阵]
int GetNumEquations()const;//获取方程个数*/




//
//线性方程组求解算法
//


//全选主元高斯消去法
BOOL CLEquations::GetRootsetGauss(CMatrix& mtxResult)
{
	int *pnJs,l,k,i,j,nIs,p,q;
	double d,t;
	
	//方程组的属性，将常数矩阵赋给解矩阵
	mtxResult=m_mtxConst;
	double *pDataCoef=m_mtxCoef.GetData();
	double *pDataConst=mtxResult.GetData();
	int n=GetNumUnkowns();
	
	//临时缓冲区，存放列数
	pnJs=new int[n];
	
	//消元
	l=1;
	for(k=0;k<=n-2;k++)
	{
		d=0.0;
		for(i=k;i<=n-1;i++)
		{
			for(j=k;j<=n-1;j++)
			{
				t=fabs(pDataCoef[i*n+j]);
				if (t>d) 
				{
					d=t;
					pnJs[k]=j;
					nIs=i;
				}
			}
			
		}
	
	if (d==0.0) l=0;
	else
	{
		if (pnJs[k]!=k)
		{
			for(i=0;i<=n-1;i++)
			{
				p=i*n+k;
				q=i*n+pnJs[k];
				t=pDataCoef[p];
				pDataCoef[p]=pDataCoef[q];
				pDataCoef[q]=t;
			}
		}
		if (nIs!=k) 
		{
			for(j=k;j<=n-1;j++)
			{
				p=k*n+j;
				q=nIs*n+j;
				t=pDataCoef[p];
				pDataCoef[p]=pDataCoef[q];
				pDataCoef[q]=t;
			}
			t=pDataCoef[k];
			pDataCoef[k]=pDataCoef[nIs];
			pDataCoef[nIs]=t;
		}
	}
	//求解失败
	if (l==0)
	{
		delete[] pnJs;
		return FALSE;
	}
	d=pDataCoef[k*n+k];
	for(j=k+1;j<=n-1;j++)
	{
		p=k*n+j;
		pDataCoef[p]=pDataCoef[p]/d;
	}
	pDataCoef[k]=pDataCoef[k]/d;
	for(i=k+1;i<=n-1;i++)
	{
		for(j=k+1;j<=n-1;j++)
		{
			p=i*n+j;
			pDataCoef[p]=pDataCoef[p]-pDataCoef[i*n+k]*pDataCoef[k*n+j];
		}
		pDataConst[i]=pDataConst[i]-pDataCoef[i*n+k]*pDataConst[k];
	}
}
  //求解失败
  d=pDataCoef[(n-1)*n+n-1];
  if (d==0.0)
{
	delete[] pnJs;
	return FALSE;
}

//求解
pDataConst[n-1]=pDataConst[n-1]/d;
for(i=n-2;i>=0;i--)
{
	t=0.0;
	for(j=i+1;j<=n-1;j++)
	{
		t=t+pDataCoef[i*n+j]*pDataConst[j];
	}
	pDataConst[i]=pDataConst[i]-t;
}
//调整解的位置
pnJs[n-1]=n-1;
for(k=n-1;k>=0;k--)
{
	if (pnJs[k]!=k)
	{
		t=pDataConst[k];
		pDataConst[k]=pDataConst[pnJs[k]];
		pDataConst[pnJs[k]]=t;
	}
}

//清理内存
delete [] pnJs;


return TRUE;

}