nlequation.cpp

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

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

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

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

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

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

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

//////////////////////////////////////////////////////////////////////
// 求非线性方程实根的对分法
//
// 调用时，须覆盖计算方程左端函数f(x)值的虚函数double Func(double x)
//
// 参数：
// 1. int nNumRoots - 在[xStart, xEnd]内实根个数的预估值
// 2. double x[] - 一维数组，长度为m。返回在区间[xStart, xEnd]内搜索到
//                 的实根，实根个数由函数值返回
// 3. double xStart - 求根区间的左端点
// 4. double xEnd - 求根区间的右端点
// 5. double dblStep - 搜索求根时采用的步长
// 6. double eps - 精度控制参数，默认值为0.000001
//
// 返回值：int 型，求得的实根的数目
//////////////////////////////////////////////////////////////////////
int CNLequation::GetRootBisect(int nNumRoots, double x[], double xStart, double xEnd, double dblStep, double eps /*= 0.000001*/)
{
	int n,js;
    double z,y,z1,y1,z0,y0;

	// 根的个数清0
    n = 0; 

	// 从左端点开始搜索
	z = xStart; 
	y = Func(z);

	// 循环求解
    while ((z<=xEnd+dblStep/2.0)&&(n!=nNumRoots))
    { 
		if (fabs(y)<eps)
        { 
			n=n+1; 
			x[n-1]=z;
            z=z+dblStep/2.0; 
			y=Func(z);
        }
        else
        { 
			z1=z+dblStep; 
			y1=Func(z1);
            
			if (fabs(y1)<eps)
            { 
				n=n+1; 
				x[n-1]=z1;
                z=z1+dblStep/2.0; 
				y=Func(z);
            }
            else if (y*y1>0.0)
            { 
				y=y1; 
				z=z1;
			}
            else
            { 
				js=0;
                while (js==0)
                { 
					if (fabs(z1-z)<eps)
                    { 
						n=n+1; 
						x[n-1]=(z1+z)/2.0;
                        z=z1+dblStep/2.0; y=Func(z);
                        js=1;
                    }
                    else
                    { 
						z0=(z1+z)/2.0; 
						y0=Func(z0);
                        if (fabs(y0)<eps)
                        { 
							x[n]=z0; 
							n=n+1; 
							js=1;
                            z=z0+dblStep/2.0; 
							y=Func(z);
                        }
                        else if ((y*y0)<0.0)
                        { 
							z1=z0; 
							y1=y0;
						}
                        else 
						{ 
							z=z0; 
							y=y0;
						}
                    }
                }
            }
        }
    }
    
	// 返回实根的数目
	return(n);
}

//////////////////////////////////////////////////////////////////////
// 求非线性方程一个实根的牛顿法
//
// 调用时，须覆盖计算方程左端函数f(x)及其一阶导数f'(x)值的虚函数
//                  void Func(double x, double y[])
//         y(0) 返回f(x)的值
//         y(1) 返回f'(x)的值
//
// 参数：
// 1. double *x - 传入迭代初值（猜测解），返回在区间求得的一个实根
// 2. int nMaxIt - 递归次数，默认值为60
// 3. double eps - 精度控制参数，默认值为0.000001
//
// 返回值：BOOL 型，求解是否成功
//////////////////////////////////////////////////////////////////////
BOOL CNLequation::GetRootNewton(double* x, int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
{ 
    int l;
    double y[2],d,p,x0,x1;

	// 条件值
    l=nMaxIt; 
	x0=*x;
    Func(x0,y);
    
	// 求解，控制精度
	d=eps+1.0;
    while ((d>=eps)&&(l!=0))
    { 
		if (y[1] == 0.0)
			return FALSE;

        x1=x0-y[0]/y[1];
        Func(x1,y);
        
		d=fabs(x1-x0); 
		p=fabs(y[0]);
        if (p>d) 
			d=p;
        x0=x1; 
		l=l-1;
    }
    
	*x=x1;

	return TRUE;
}

//////////////////////////////////////////////////////////////////////
// 求非线性方程一个实根的埃特金迭代法
//
// 调用时，须覆盖计算方程左端函数f(x)值的虚函数double Func(double x)
//
// 参数：
// 1. double *x - 传入迭代初值（猜测解），返回在区间求得的一个实根
// 2. int nMaxIt - 递归次数，默认值为60
// 3. double eps - 精度控制参数，默认值为0.000001
//
// 返回值：BOOL 型，求解是否成功
//////////////////////////////////////////////////////////////////////
BOOL CNLequation::GetRootAitken(double* x, int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
{ 
	int flag,l;
    double u,v,x0;

	// 求解条件
    l=0; 
	x0=*x; 
	flag=0;

	// 迭代求解
    while ((flag==0)&&(l!=nMaxIt))
    { 
		l=l+1; 
        u=Func(x0); 
		v=Func(u);
        if (fabs(u-v)<eps) 
		{ 
			x0=v; 
			flag=1; 
		}
        else 
			x0=v-(v-u)*(v-u)/(v-2.0*u+x0);
    }
    
	*x=x0; 
    
	// 是否在指定的迭代次数内达到求解精度
	return (nMaxIt > l);
}

//////////////////////////////////////////////////////////////////////
// 求非线性方程一个实根的连分式解法
//
// 调用时，须覆盖计算方程左端函数f(x)值的虚函数double Func(double x)
//
// 参数：
// 1. double *x - 传入迭代初值（猜测解），返回在区间求得的一个实根
// 2. double eps - 精度控制参数，默认值为0.000001
//
// 返回值：BOOL 型，求解是否成功
//////////////////////////////////////////////////////////////////////
BOOL CNLequation::GetRootPq(double* x, double eps /*= 0.000001*/)
{ 
    int i,j,m,it,l;
    double a[10],y[10],z,h,x0,q;

	// 求解条件
    l=10; 
	q=1.0e+35; 
	x0=*x; 
	h=0.0;
    
	// 连分式求解
	while (l!=0)
    { 
		l=l-1; 
		j=0; 
		it=l;
        while (j<=7)
        { 
			if (j<=2) 
				z=x0+0.1*j;
            else 
				z=h;
			
			y[j]=Func(z);
            h=z;
            if (j==0) 
				a[0]=z;
            else
            { 
				m=0; 
				i=0;
                while ((m==0)&&(i<=j-1))
                { 
					if (fabs(h-a[i])+1.0==1.0) 
						m=1;
                    else 
						h=(y[j]-y[i])/(h-a[i]);
                     
					i=i+1;
                }
                a[j]=h;
                if (m!=0) 
					a[j]=q;
                h=0.0;
                for (i=j-1; i>=0; i--)
                { 
					if (fabs(a[i+1]+h)+1.0==1.0) 
						h=q;
                    else 
						h=-y[i]/(a[i+1]+h);
                }
                 
				h=h+a[0];
            }
             
			if (fabs(y[j])>=eps) 
				j=j+1;
            else 
			{ 
				j=10; 
				l=0;
			}
		}
        
		x0=h;
	}

    *x=h;
    
	// 是否在10阶连分式内求的实根？
	return (10>it);
}

//////////////////////////////////////////////////////////////////////
// 求实系数代数方程全部根的QR方法
//
// 参数：
// 1. int n - 多项式方程的次数
// 2. double dblCoef[] - 一维数组，长度为n+1，按降幂次序依次存放n次多项式方程的n+1个系数
// 3. double xr[] - 一维数组，长度为n，返回n个根的实部
// 4. double xi[] - 一维数组，长度为n，返回n个根的虚部
// 5. int nMaxIt - 迭代次数，默认值为 60
// 6. double eps - 精度控制参数，默认值为 0.000001
//
// 返回值：BOOL 型，求解是否成功
//////////////////////////////////////////////////////////////////////
BOOL CNLequation::GetRootQr(int n, double dblCoef[], double xr[], double xi[], int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
{ 
	// 初始化矩阵
	CMatrix mtxQ;
	mtxQ.Init(n, n);
    double *q = mtxQ.GetData();

	// 构造赫申伯格矩阵
    for (int j=0; j<=n-1; j++)
		q[j]=-dblCoef[n-j-1]/dblCoef[n];

    for (j=n; j<=n*n-1; j++)
		q[j]=0.0;

    for (int i=0; i<=n-2; i++)
		q[(i+1)*n+i]=1.0;

	// 求赫申伯格矩阵的特征值和特征向量，即为方程的解
    if (mtxQ.HBergEigenv(xr, xi, nMaxIt, eps))
		return TRUE;

	return FALSE;
}

//////////////////////////////////////////////////////////////////////
// 求实系数代数方程全部根的牛顿下山法
//
// 参数：
// 1. int n - 多项式方程的次数
// 2. double dblCoef[] - 一维数组，长度为n+1，按降幂次序依次存放n次多项式方程的n+1个系数
// 3. double xr[] - 一维数组，长度为n，返回n个根的实部
// 4. double xi[] - 一维数组，长度为n，返回n个根的虚部
//
// 返回值：BOOL 型，求解是否成功
//////////////////////////////////////////////////////////////////////
BOOL CNLequation::GetRootNewtonDownHill(int n, double dblCoef[], double xr[], double xi[])
{ 
	int m,i,jt,k,is,it;
    double t,x,y,x1,y1,dx,dy,p,q,w,dd,dc,c;
    double g,u,v,pq,g1,u1,v1;
    
	// 初始判断
    m=n;
    while ((m>0)&&(fabs(dblCoef[m])+1.0==1.0)) 
		m=m-1;

	// 求解失败
    if (m<=0)
		return FALSE;

    for (i=0; i<=m; i++)
		dblCoef[i]=dblCoef[i]/dblCoef[m];
    
	for (i=0; i<=m/2; i++)
    { 
		w=dblCoef[i]; 
		dblCoef[i]=dblCoef[m-i]; 
		dblCoef[m-i]=w;
	}
    
	// 迭代求解
	k=m; 
	is=0; 
	w=1.0;
    jt=1;
    while (jt==1)
    { 
		pq=fabs(dblCoef[k]);
		while (pq<1.0e-12)
        { 
			xr[k-1]=0.0; 
			xi[k-1]=0.0; 
			k=k-1;
            if (k==1)
            { 
				xr[0]=-dblCoef[1]*w/dblCoef[0]; 
				xi[0]=0.0;
                
				return TRUE;
            }
            
			pq=fabs(dblCoef[k]);
        }
	
		q=log(pq); 
		q=q/(1.0*k); 
		q=exp(q);
        p=q; 
		w=w*p;
        for (i=1; i<=k; i++)
        { 
			dblCoef[i]=dblCoef[i]/q; 
			q=q*p;
		}
        
		x=0.0001; 
		x1=x; 
		y=0.2; 
		y1=y; 
		dx=1.0;
        g=1.0e+37; 
l40:
        u=dblCoef[0]; v=0.0;
        for (i=1; i<=k; i++)
        { 
			p=u*x1; 
			q=v*y1;
            pq=(u+v)*(x1+y1);
            u=p-q+dblCoef[i]; 
			v=pq-p-q;
        }
        
		g1=u*u+v*v;
        if (g1>=g)
        { 
			if (is!=0)
            { 
				it=1;
                g65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,&k,&is,&it);
                if (it==0) 
					goto l40;
            }
            else
            { 
				g60(&t,&x,&y,&x1,&y1,&dx,&dy,&p,&q,&k,&it);
                if (t>=1.0e-03) 
					goto l40;
                
				if (g>1.0e-18)
                { 
					it=0;
                    g65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,&k,&is,&it);
                    if (it==0) 
						goto l40;
                }
            }
            
			g90(xr,xi,dblCoef,&x,&y,&p,&q,&w,&k);
        }
        else
        { 
			g=g1; 
			x=x1; 
			y=y1; 
			is=0;
            if (g<=1.0e-22)
				g90(xr,xi,dblCoef,&x,&y,&p,&q,&w,&k);
            else
            { 
				u1=k*dblCoef[0]; 
				v1=0.0;
                for (i=2; i<=k; i++)
                {