nlequation.cpp

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

			}
        }
    }
    
	*x=xx;
}

//////////////////////////////////////////////////////////////////////
// 内部函数
//////////////////////////////////////////////////////////////////////
double CNLequation::rnd(double* r)
{
	int m;
    double s,u,v,p;
    
	s=65536.0; 
	u=2053.0; 
	v=13849.0;
    m=(int)(*r/s); 
	*r=*r-m*s;
    *r=u*(*r)+v; 
	m=(int)(*r/s);
    *r=*r-m*s; p=*r/s;
    
	return(p);
}

//////////////////////////////////////////////////////////////////////
// 求实函数或复函数方程一个复根的蒙特卡洛法
//
// 调用时，须覆盖计算方程左端函数的模值||f(x, y)||的虚函数
//           double Func(double x, double y)
//
// 参数：
// 1. double* x - 传入初值（猜测解）的实部，返回求得的根的实部
// 2. double* y - 传入初值（猜测解）的虚部，返回求得的根的虚部
// 3. double xStart - 均匀分布的端点初值
// 4. int nControlB - 控制参数
// 5. double eps - 控制精度，默认值为0.000001
// 6. double xi[] - 一维数组，长度为n，返回n个根的虚部
//
// 返回值：无
//////////////////////////////////////////////////////////////////////
void CNLequation::GetRootMonteCarlo(double* x, double* y, double xStart, int nControlB, double eps /*= 0.000001*/)
{ 
    int k;
    double xx,yy,a,r,z,x1,y1,z1;

	// 求解条件与初值
    a=xStart; 
	k=1; 
	r=1.0; 
	xx=*x; 
	yy=*y;
    z=Func(xx,yy);
    
	// 精度控制求解
	while (a>=eps)
    { 
		x1=-a+2.0*a*rnd(&r); 
		x1=xx+x1; 
        y1=-a+2.0*a*rnd(&r); 
		y1=yy+y1;
        
		z1=Func(x1,y1);
        
		k=k+1;
        if (z1>=z)
        { 
			if (k>nControlB) 
			{ 
				k=1; 
				a=a/2.0; 
			}
		}
        else
        { 
			k=1; 
			xx=x1; 
			yy=y1;  
			z=z1;
            if (z<eps)
            { 
				*x=xx; 
				*y=yy; 
				return; 
			}
        }
    }
    
	*x=xx; 
	*y=yy; 
}

//////////////////////////////////////////////////////////////////////
// 求非线性方程组一组实根的梯度法
//
// 调用时，须覆盖计算方程左端函数f(x)值及其偏导数值的虚函数
//         double Func(double x[], double y[])
//
// 参数：
// 1. int n - 方程的个数，也是未知数的个数
// 2. double x[] - 一维数组，长度为n，存放一组初值x0, x1, …, xn-1，
//                 返回时存放方程组的一组实根
// 3. int nMaxIt - 迭代次数，默认值为60
// 4. double eps - 控制精度，默认值为0.000001
//
// 返回值：BOOL 型，求解是否成功
//////////////////////////////////////////////////////////////////////
BOOL CNLequation::GetRootsetGrad(int n, double x[], int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
{ 
    int l,j;
    double f,d,s,*y;

    y = new double[n];

    l=nMaxIt;
    f=Func(x,y);

	// 控制精度，迭代求解
    while (f>=eps)
    { 
		l=l-1;
        if (l==0) 
		{ 
			delete[] y; 
			return TRUE;
		}
        
		d=0.0;
        for (j=0; j<=n-1; j++) 
			d=d+y[j]*y[j];
        if (d+1.0==1.0) 
		{ 
			delete[] y; 
			return FALSE;
		}
        
		s=f/d;
        for (j=0; j<=n-1; j++) 
			x[j]=x[j]-s*y[j];
        
		f=Func(x,y);
    }
    
	delete[] y; 

	// 是否在有效迭代次数内达到精度
	return (nMaxIt>l);
}

//////////////////////////////////////////////////////////////////////
// 求非线性方程组一组实根的拟牛顿法
//
// 调用时，须覆盖计算方程左端函数f(x)值及其偏导数值的虚函数
//         double Func(double x[], double y[])
//
// 参数：
// 1. int n - 方程的个数，也是未知数的个数
// 2. double x[] - 一维数组，长度为n，存放一组初值x0, x1, …, xn-1，
//                 返回时存放方程组的一组实根
// 3. double t - 控制h大小的变量，0<t<1
// 4. double h - 增量初值
// 3. int nMaxIt - 迭代次数，默认值为60
// 4. double eps - 控制精度，默认值为0.000001
//
// 返回值：BOOL 型，求解是否成功
//////////////////////////////////////////////////////////////////////
BOOL CNLequation::GetRootsetNewton(int n, double x[], double t, double h, int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
{ 
    int i,j,l;
    double am,z,beta,d,*y,*a,*b;

    y = new double[n];
	
	// 构造矩阵
	CMatrix mtxCoef(n, n);
	CMatrix mtxConst(n, 1);
    a = mtxCoef.GetData();
    b = mtxConst.GetData();

	// 迭代求解
    l=nMaxIt; 
	am=1.0+eps;
    while (am>=eps)
    { 
		Func(x,b);
        
		am=0.0;
        for (i=0; i<=n-1; i++)
        { 
			z=fabs(b[i]);
            if (z>am) 
				am=z;
        }
        
		if (am>=eps)
        { 
			l=l-1;
            if (l==0)
            { 
				delete[] y;
                return FALSE;
            }
            
			for (j=0; j<=n-1; j++)
            { 
				z=x[j]; 
				x[j]=x[j]+h;
                
				Func(x,y);
                
				for (i=0; i<=n-1; i++) 
					a[i*n+j]=y[i];
                
				x[j]=z;
            }
            
			// 调用全选主元高斯消元法
			CLEquations leqs(mtxCoef, mtxConst);
			CMatrix mtxResult;
			if (! leqs.GetRootsetGauss(mtxResult))
			{
				delete[] y;
				return FALSE;
			}

			mtxConst = mtxResult;
			b = mtxConst.GetData();

			beta=1.0;
            for (i=0; i<=n-1; i++) 
				beta=beta-b[i];
            
			if (beta == 0.0)
            { 
				delete[] y;
                return FALSE;
            }
            
			d=h/beta;
            for (i=0; i<=n-1; i++) 
				x[i]=x[i]-d*b[i];
            
			h=t*h;
        }
    }
    
	delete[] y;
    
	// 是否在有效迭代次数内达到精度
	return (nMaxIt>l);
}

//////////////////////////////////////////////////////////////////////
// 求非线性方程组最小二乘解的广义逆法
//
// 调用时，1. 须覆盖计算方程左端函数f(x)值及其偏导数值的虚函数
//             double Func(double x[], double y[])
//         2. 须覆盖计算雅可比矩阵函数的虚函数
//             double FuncMJ(int n, double x[], double y[])
//
// 参数：
// 1. int m - 方程的个数
// 2. int n - 未知数的个数
// 3. double x[] - 一维数组，长度为n，存放一组初值x0, x1, …, xn-1，要求
//                 不全为0，返回时存放方程组的最小二乘解，当m=n时，即是
//                 非线性方程组的解
// 4. double eps1 - 最小二乘解的精度控制精度，默认值为0.000001
// 5. double eps2 - 奇异值分解的精度控制精度，默认值为0.000001
//
// 返回值：BOOL 型，求解是否成功
//////////////////////////////////////////////////////////////////////
BOOL CNLequation::GetRootsetGinv(int m, int n, double x[], double eps1 /*= 0.000001*/, double eps2 /*= 0.000001*/)
{ 
	int i,j,k,l,kk,jt;
    double y[10],b[10],alpha,z,h2,y1,y2,y3,y0,h1;
    double *p,*d,*dx,*w;
    
	// 控制参数
	int ka = max(m, n)+1;
    w = new double[ka];
    
	// 设定迭代次数为60，迭代求解
	l=60; 
	alpha=1.0;
    while (l>0)
    { 
		CMatrix mtxP(m, n), mtxD(m, 1);
		p = mtxP.GetData();
		d = mtxD.GetData();

		Func(x,d);
        FuncMJ(n,x,p);

		// 构造线性方程组
		CLEquations leqs(mtxP, mtxD);
		// 临时矩阵
		CMatrix mtxAP, mtxU, mtxV;
		// 解矩阵
		CMatrix mtxDX;
		// 基于广义逆的最小二乘解
		if (! leqs.GetRootsetGinv(mtxDX, mtxAP, mtxU, mtxV, eps2))
        { 
			delete[] w;			
			return FALSE;
        }
        
		dx = mtxDX.GetData();

		j=0; 
		jt=1; 
		h2=0.0;
        while (jt==1)
        { 
			jt=0;
            if (j<=2) 
				z=alpha+0.01*j;
            else 
				z=h2;
            
			for (i=0; i<=n-1; i++) 
				w[i]=x[i]-z*dx[i];
            
			Func(w,d);
            
			y1=0.0;
            for (i=0; i<=m-1; i++) 
				y1=y1+d[i]*d[i];
            for (i=0; i<=n-1; i++)
				w[i]=x[i]-(z+0.00001)*dx[i];
            
			Func(w,d);
            
			y2=0.0;
            for (i=0; i<=m-1; i++) 
				y2=y2+d[i]*d[i];
            
			y0=(y2-y1)/0.00001;
            
			if (fabs(y0)>1.0e-10)
            { 
				h1=y0; h2=z;
                if (j==0) 
				{ 
					y[0]=h1; 
					b[0]=h2;
				}
                else
                { 
					y[j]=h1; 
					kk=0; 
					k=0;
                    while ((kk==0)&&(k<=j-1))
                    { 
						y3=h2-b[k];
                        if (fabs(y3)+1.0==1.0) 
							kk=1;
                        else 
							h2=(h1-y[k])/y3;
                        
						k=k+1;
                    }
                    
					b[j]=h2;
                    if (kk!=0) 
						b[j]=1.0e+35;
                    
					h2=0.0;
                    for (k=j-1; k>=0; k--)
						h2=-y[k]/(b[k+1]+h2);
                    
					h2=h2+b[0];
                }
                
				j=j+1;
                if (j<=7) 
					jt=1;
                else 
					z=h2;
            }
        }
        
		alpha=z; 
		y1=0.0; 
		y2=0.0;
        for (i=0; i<=n-1; i++)
        { 
			dx[i]=-alpha*dx[i];
            x[i]=x[i]+dx[i];
            y1=y1+fabs(dx[i]);
            y2=y2+fabs(x[i]);
        }
        
		// 求解成功
		if (y1<eps1*y2)
        { 
			delete[] w;			
			return TRUE;
        }
        
		l=l-1;
    }
    
	delete[] w;			
	
	// 求解失败
	return FALSE;
}

//////////////////////////////////////////////////////////////////////
// 求非线性方程组一组实根的蒙特卡洛法
//
// 调用时，须覆盖计算方程左端模函数值||F||的虚函数
//         double Func(int n, double x[])
//         其返回值为Sqr(f1*f1 + f2*f2 + … + fn*fn)
//
// 参数：
// 1. double x[] - 一维数组，长度为n，存放一组初值，返回时存放方程的一组实根
// 2. double xStart - 均匀分布的端点初值
// 3. int nControlB - 控制参数
// 4. double eps - 控制精度，默认值为0.000001
//
// 返回值：无
//////////////////////////////////////////////////////////////////////
void CNLequation::GetRootsetMonteCarlo(int n, double x[], double xStart, int nControlB, double eps /*= 0.000001*/)
{ 
	int k,i;
    double a,r,*y,z,z1;

    y = new double[n];
    
	// 初值
	a=xStart; 
	k=1; 
	r=1.0; 

	z = Func(n, x);

	// 用精度控制迭代求解
    while (a>=eps)
    { 
		for (i=0; i<=n-1; i++)
			y[i]=-a+2.0*a*rnd(&r)+x[i];
        
		z1=Func(n,y);
        
		k=k+1;
        if (z1>=z)
        { 
			if (k>nControlB) 
			{ 
				k=1; 
				a=a/2.0; 
			}
		}
        else
        { 
			k=1; 
            for (i=0; i<=n-1; i++) 
				x[i]=y[i];
            
			// 求解成功
			z=z1;
            if (z<eps)  
			{
				delete[] y;
				return;
			}
        }
    }

	delete[] y;
}