lequations.cpp

科学与工程数值算法求解方程组的类

//////////////////////////////////////////////////////////////////////
// LEquations.cpp
//
// 求解线性方程组的类 CLEquations 的实现代码
//
// 周长发编制, 2002/8
//////////////////////////////////////////////////////////////////////

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

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

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

//////////////////////////////////////////////////////////////////////
// 基本构造函数
//////////////////////////////////////////////////////////////////////
CLEquations::CLEquations()
{
}

//////////////////////////////////////////////////////////////////////
// 指定系数和常数构造函数
//
// 参数：
// 1. const CMatrix& mtxCoef - 指定的系数矩阵
// 2. const CMatrix& mtxConst - 指定的常数矩阵
//////////////////////////////////////////////////////////////////////
CLEquations::CLEquations(const CMatrix& mtxCoef, const CMatrix& mtxConst)
{
	ASSERT(Init(mtxCoef, mtxConst));
}

//////////////////////////////////////////////////////////////////////
// 析构函数
//////////////////////////////////////////////////////////////////////
CLEquations::~CLEquations()
{
}

//////////////////////////////////////////////////////////////////////
// 初始化函数
//
// 参数：
// 1. const CMatrix& mtxCoef - 指定的系数矩阵
// 2. const CMatrix& mtxConst - 指定的常数矩阵
//
// 返回值：BOOL 型，初始化是否成功
//////////////////////////////////////////////////////////////////////
BOOL CLEquations::Init(const CMatrix& mtxCoef, const CMatrix& mtxConst)
{
	if (mtxCoef.GetNumRows() != mtxConst.GetNumRows())
		return FALSE;

	m_mtxCoef = mtxCoef;
	m_mtxConst = mtxConst;

	return TRUE;
}

//////////////////////////////////////////////////////////////////////
// 获取系数矩阵
//
// 参数：无
//
// 返回值：CMatrix 型，返回系数矩阵
//////////////////////////////////////////////////////////////////////
inline CMatrix CLEquations::GetCoefMatrix() const
{
	return m_mtxCoef;
}

//////////////////////////////////////////////////////////////////////
// 获取常数矩阵
//
// 参数：无
//
// 返回值：CMatrix 型，返回常数矩阵
//////////////////////////////////////////////////////////////////////
inline CMatrix CLEquations::GetConstMatrix() const
{
	return m_mtxConst;
}

//////////////////////////////////////////////////////////////////////
// 获取方程个数
//
// 参数：无
//
// 返回值：int 型，返回方程组方程的个数
//////////////////////////////////////////////////////////////////////
inline int	CLEquations::GetNumEquations() const
{
	return GetCoefMatrix().GetNumRows();
}

//////////////////////////////////////////////////////////////////////
// 获取未知数个数
//
// 参数：无
//
// 返回值：int 型，返回方程组未知数的个数
//////////////////////////////////////////////////////////////////////
inline int	CLEquations::GetNumUnknowns() const
{
	return GetCoefMatrix().GetNumColumns();
}

//////////////////////////////////////////////////////////////////////
// 全选主元高斯消去法
//
// 参数：
// 1. CMatrix& mtxResult - CMatrix引用对象，返回方程组的解
//
// 返回值：BOOL 型，方程组求解是否成功
//////////////////////////////////////////////////////////////////////
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 = GetNumUnknowns();

	// 临时缓冲区，存放列数
    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=pDataConst[k]; 
				pDataConst[k]=pDataConst[nIs]; 
				pDataConst[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;
		}
        
		pDataConst[k]=pDataConst[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;
}

//////////////////////////////////////////////////////////////////////
// 全选主元高斯－约当消去法
//
// 参数：
// 1. CMatrix& mtxResult - CMatrix引用对象，返回方程组的解
//
// 返回值：BOOL 型，方程组求解是否成功
//////////////////////////////////////////////////////////////////////
BOOL CLEquations::GetRootsetGaussJordan(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 = GetNumUnknowns();
	int m = m_mtxConst.GetNumColumns();

	// 临时缓冲区，存放变换的列数
    pnJs = new int[n];

	// 消元
    l=1;
    for (k=0;k<=n-1;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+1.0==1.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;
                }
                
				for (j=0;j<=m-1;j++)
                { 
					p=k*m+j; 
					q=nIs*m+j;
                    t=pDataConst[p]; 
					pDataConst[p]=pDataConst[q]; 
					pDataConst[q]=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;
		}
        
		for (j=0;j<=m-1;j++)
        { 
			p=k*m+j; 
			pDataConst[p]=pDataConst[p]/d;
		}
        
		for (j=k+1;j<=n-1;j++)
        {
			for (i=0;i<=n-1;i++)
            { 
				p=i*n+j;
				if (i!=k)
					pDataCoef[p]=pDataCoef[p]-pDataCoef[i*n+k]*pDataCoef[k*n+j];
            }
		}

        for (j=0;j<=m-1;j++)
        {
			for (i=0;i<=n-1;i++)
			{ 
				p=i*m+j;
				if (i!=k)
				pDataConst[p]=pDataConst[p]-pDataCoef[i*n+k]*pDataConst[k*m+j];
			}
		}
    }
    
	// 调整
	for (k=n-1;k>=0;k--)
    {
		if (pnJs[k]!=k)
		{
			for (j=0;j<=m-1;j++)
			{ 
				p=k*m+j; 
				q=pnJs[k]*m+j;
				t=pDataConst[p]; 
				pDataConst[p]=pDataConst[q]; 
				pDataConst[q]=t;
			}
		}
	}

	// 清理内存
    delete[] pnJs;
    
	return TRUE;
}

//////////////////////////////////////////////////////////////////////
// 复系数方程组的全选主元高斯消去法
//
// 参数：
// 1. const CMatrix& mtxCoefImag - 系数矩阵的虚部矩阵
// 2. const CMatrix& mtxConstImag - 常数矩阵的虚部矩阵
// 3. CMatrix& mtxResult - CMatrix引用对象，返回方程组解矩阵的实部矩阵
// 4. CMatrix& mtxResultImag - CMatrix引用对象，返回方程组解矩阵的虚部矩阵
//
// 返回值：BOOL 型，方程组求解是否成功
//////////////////////////////////////////////////////////////////////
BOOL CLEquations::GetRootsetGauss(const CMatrix& mtxCoefImag, const CMatrix& mtxConstImag, CMatrix& mtxResult, CMatrix& mtxResultImag)
{ 
	int *pnJs,l,k,i,j,nIs,u,v;
    double p,q,s,d;

	// 方程组的属性，将常数矩阵赋给解矩阵
	mtxResult = m_mtxConst;
	mtxResultImag = mtxConstImag;
	double *pDataCoef = m_mtxCoef.GetData();
	double *pDataConst = mtxResult.GetData();
	double *pDataCoefImag = mtxCoefImag.GetData();
	double *pDataConstImag = mtxResultImag.GetData();
	int n = GetNumUnknowns();
	int m = m_mtxConst.GetNumColumns();

	// 临时缓冲区，存放变换的列数
    pnJs = new int[n];
    
	// 消元
	for (k=0;k<=n-2;k++)
    { 
		d=0.0;
        for (i=k;i<=n-1;i++)
        {
			for (j=k;j<=n-1;j++)
			{ 
				u=i*n+j;
				p=pDataCoef[u]*pDataCoef[u]+pDataCoefImag[u]*pDataCoefImag[u];
				if (p>d) 
				{
					d=p;
					pnJs[k]=j;
					nIs=i;
				}
			}
        }
        
		// 求解失败
		if (d == 0.0)
		{
			delete[] pnJs;
			return FALSE;
        }
        
		if (nIs!=k)
        { 
			for (j=k;j<=n-1;j++)
            { 
				u=k*n+j; 
				v=nIs*n+j;
                p=pDataCoef[u]; 
				pDataCoef[u]=pDataCoef[v]; 
				pDataCoef[v]=p;
                p=pDataCoefImag[u]; 
				pDataCoefImag[u]=pDataCoefImag[v]; 
				pDataCoefImag[v]=p;
            }
            
			p=pDataConst[k]; 
			pDataConst[k]=pDataConst[nIs]; 
			pDataConst[nIs]=p;
            p=pDataConstImag[k]; 
			pDataConstImag[k]=pDataConstImag[nIs]; 
			pDataConstImag[nIs]=p;
        }
        
		if (pnJs[k]!=k)
        {
			for (i=0;i<=n-1;i++)
            { 
				u=i*n+k; 
				v=i*n+pnJs[k];
				p=pDataCoef[u]; 
				pDataCoef[u]=pDataCoef[v]; 
				pDataCoef[v]=p;
				p=pDataCoefImag[u]; 
				pDataCoefImag[u]=pDataCoefImag[v]; 
				pDataCoefImag[v]=p;
            }
		}

        v=k*n+k;
        for (j=k+1;j<=n-1;j++)
        { 
			u=k*n+j;
            p=pDataCoef[u]*pDataCoef[v]; 
			q=-pDataCoefImag[u]*pDataCoefImag[v];
            s=(pDataCoef[v]-pDataCoefImag[v])*(pDataCoef[u]+pDataCoefImag[u]);
            pDataCoef[u]=(p-q)/d; 
			pDataCoefImag[u]=(s-p-q)/d;
        }
        
		p=pDataConst[k]*pDataCoef[v]; 
		q=-pDataConstImag[k]*pDataCoefImag[v];
        s=(pDataCoef[v]-pDataCoefImag[v])*(pDataConst[k]+pDataConstImag[k]);
        pDataConst[k]=(p-q)/d; 
		pDataConstImag[k]=(s-p-q)/d;

        for (i=k+1;i<=n-1;i++)
        { 
			u=i*n+k;
            for (j=k+1;j<=n-1;j++)
            { 
				v=k*n+j; 
				l=i*n+j;
                p=pDataCoef[u]*pDataCoef[v]; 
				q=pDataCoefImag[u]*pDataCoefImag[v];
                s=(pDataCoef[u]+pDataCoefImag[u])*(pDataCoef[v]+pDataCoefImag[v]);
                pDataCoef[l]=pDataCoef[l]-p+q;
                pDataCoefImag[l]=pDataCoefImag[l]-s+p+q;
            }
            
			p=pDataCoef[u]*pDataConst[k]; 
			q=pDataCoefImag[u]*pDataConstImag[k];
            s=(pDataCoef[u]+pDataCoefImag[u])*(pDataConst[k]+pDataConstImag[k]);
            pDataConst[i]=pDataConst[i]-p+q; 
			pDataConstImag[i]=pDataConstImag[i]-s+p+q;
        }
    }
    
	u=(n-1)*n+n-1;
    d=pDataCoef[u]*pDataCoef[u]+pDataCoefImag[u]*pDataCoefImag[u];

	// 求解失败
    if (d == 0.0)
	{
		delete[] pnJs;
		return FALSE;
    }

	// 求解
    p=pDataCoef[u]*pDataConst[n-1]; q=-pDataCoefImag[u]*pDataConstImag[n-1];
    s=(pDataCoef[u]-pDataCoefImag[u])*(pDataConst[n-1]+pDataConstImag[n-1]);
    pDataConst[n-1]=(p-q)/d; pDataConstImag[n-1]=(s-p-q)/d;

    for (i=n-2;i>=0;i--)
    {
		for (j=i+1;j<=n-1;j++)
		{ 
			u=i*n+j;
			p=pDataCoef[u]*pDataConst[j]; 
			q=pDataCoefImag[u]*pDataConstImag[j];
			s=(pDataCoef[u]+pDataCoefImag[u])*(pDataConst[j]+pDataConstImag[j]);
			pDataConst[i]=pDataConst[i]-p+q;
			pDataConstImag[i]=pDataConstImag[i]-s+p+q;
		}
    }

	// 调整位置
    pnJs[n-1]=n-1;
    for (k=n-1;k>=0;k--)
    {
		if (pnJs[k]!=k)
        { 
			p=pDataConst[k]; 
			pDataConst[k]=pDataConst[pnJs[k]]; 
			pDataConst[pnJs[k]]=p;
			p=pDataConstImag[k]; 
			pDataConstImag[k]=pDataConstImag[pnJs[k]]; 
			pDataConstImag[pnJs[k]]=p;
        }
	}

	// 清理内存
	delete[] pnJs;

	return TRUE;
}

//////////////////////////////////////////////////////////////////////
// 复系数方程组的全选主元高斯－约当消去法
//
// 参数：
// 1. const CMatrix& mtxCoefImag - 系数矩阵的虚部矩阵
// 2. const CMatrix& mtxConstImag - 常数矩阵的虚部矩阵
// 3. CMatrix& mtxResult - CMatrix引用对象，返回方程组解矩阵的实部矩阵