nlequation.cpp

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

					p=u1*x; 
					q=v1*y; 
					pq=(u1+v1)*(x+y);
                    u1=p-q+(k-i+1)*dblCoef[i-1];
                    v1=pq-p-q;
                }
                
				p=u1*u1+v1*v1;
                if (p<=1.0e-20)
                { 
					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
                { 
					dx=(u*u1+v*v1)/p;
                    dy=(u1*v-v1*u)/p;
                    t=1.0+4.0/k;
                    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);
                }
            }
        }
        
		if (k==1) 
			jt=0;
        else 
			jt=1;
    }
    
	return TRUE;
}


//////////////////////////////////////////////////////////////////////
// 内部函数
//////////////////////////////////////////////////////////////////////
void CNLequation::g60(double* t,double* x,double* y,double* x1,double* y1,double* dx,double* dy,double* p,double* q,int* k,int* it)
{ 
	*it=1;
    while (*it==1)
    { 
		*t=*t/1.67; 
		*it=0;
        *x1=*x-(*t)*(*dx);
        *y1=*y-(*t)*(*dy);
        if (*k>=50)
		{ 
			*p=sqrt((*x1)*(*x1)+(*y1)*(*y1));
            *q=exp(85.0/(*k));
            if (*p>=*q) 
				*it=1;
        }
    }
}

//////////////////////////////////////////////////////////////////////
// 内部函数
//////////////////////////////////////////////////////////////////////
void CNLequation::g90(double xr[],double xi[],double dblCoef[],double* x,double* y,double* p,double* q,double* w,int* k)
{ 
	int i;
    
	if (fabs(*y)<=1.0e-06)
    { 
		*p=-(*x); 
		*y=0.0; 
		*q=0.0;
	}
    else
    { 
		*p=-2.0*(*x); 
		*q=(*x)*(*x)+(*y)*(*y);
        xr[*k-1]=(*x)*(*w);
        xi[*k-1]=-(*y)*(*w);
        *k=*k-1;
    }
    
	for (i=1; i<=*k; i++)
    { 
		dblCoef[i]=dblCoef[i]-dblCoef[i-1]*(*p);
        dblCoef[i+1]=dblCoef[i+1]-dblCoef[i-1]*(*q);
    }
    
	xr[*k-1]=(*x)*(*w); 
	xi[*k-1]=(*y)*(*w);
    *k=*k-1;
    if (*k==1)
    { 
		xr[0]=-dblCoef[1]*(*w)/dblCoef[0]; 
		xi[0]=0.0;
	}
}

//////////////////////////////////////////////////////////////////////
// 内部函数
//////////////////////////////////////////////////////////////////////
void CNLequation::g65(double* x,double* y,double* x1,double* y1,double* dx,double* dy,double* dd,double* dc,double* c,int* k,int* is,int* it)
{ 
	if (*it==0)
    { 
		*is=1;
        *dd=sqrt((*dx)*(*dx)+(*dy)*(*dy));
        if (*dd>1.0) 
			*dd=1.0;
        *dc=6.28/(4.5*(*k)); 
		*c=0.0;
    }
    
	while(TRUE)
    { 
		*c=*c+(*dc);
        *dx=(*dd)*cos(*c); 
		*dy=(*dd)*sin(*c);
        *x1=*x+*dx; 
		*y1=*y+*dy;
        if (*c<=6.29)
        { 
			*it=0; 
			return;
		}
        
		*dd=*dd/1.67;
        if (*dd<=1.0e-07)
        { 
			*it=1; 
			return;
		}
        
		*c=0.0;
    }
}

//////////////////////////////////////////////////////////////////////
// 求复系数代数方程全部根的牛顿下山法
//
// 参数：
// 1. int n - 多项式方程的次数
// 2. double ar[] - 一维数组，长度为n+1，按降幂次序依次存放n次多项式方程的n+1个系数的实部
// 3. double ai[] - 一维数组，长度为n+1，按降幂次序依次存放n次多项式方程的n+1个系数的虚部
// 4. double xr[] - 一维数组，长度为n，返回n个根的实部
// 5. double xi[] - 一维数组，长度为n，返回n个根的虚部
//
// 返回值：BOOL 型，求解是否成功
//////////////////////////////////////////////////////////////////////
BOOL CNLequation::GetRootNewtonDownHill(int n, double ar[], double ai[], 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;
    p=sqrt(ar[m]*ar[m]+ai[m]*ai[m]);
    while ((m>0)&&(p+1.0==1.0))
    {  
		m=m-1;
        p=sqrt(ar[m]*ar[m]+ai[m]*ai[m]);
    }
    
	// 求解失败
	if (m<=0)
		return FALSE;

    for (i=0; i<=m; i++)
    { 
		ar[i]=ar[i]/p; 
		ai[i]=ai[i]/p;
	}
    
	for (i=0; i<=m/2; i++)
    { 
		w=ar[i]; 
		ar[i]=ar[m-i]; 
		ar[m-i]=w;
        w=ai[i]; 
		ai[i]=ai[m-i]; 
		ai[m-i]=w;
    }
    
	// 迭代求解
	k=m; 
	is=0; 
	w=1.0;
    jt=1;
    while (jt==1)
    { 
		pq=sqrt(ar[k]*ar[k]+ai[k]*ai[k]);
		while (pq<1.0e-12)
        { 
			xr[k-1]=0.0; 
			xi[k-1]=0.0; 
			k=k-1;
            if (k==1)
            { 
				p=ar[0]*ar[0]+ai[0]*ai[0];
                xr[0]=-w*(ar[0]*ar[1]+ai[0]*ai[1])/p;
                xi[0]=w*(ar[1]*ai[0]-ar[0]*ai[1])/p;
                
				return TRUE;
            }
            
			pq=sqrt(ar[k]*ar[k]+ai[k]*ai[k]);
        }
		
		q=log(pq); 
		q=q/(1.0*k); 
		q=exp(q);
        p=q; 
		w=w*p;
        for (i=1; i<=k; i++)
        { 
			ar[i]=ar[i]/q; 
			ai[i]=ai[i]/q; 
			q=q*p;
		}
        
		x=0.0001; 
		x1=x; 
		y=0.2; 
		y1=y; 
		dx=1.0;
        g=1.0e+37; 
l40:
        u=ar[0]; 
		v=ai[0];
        for (i=1; i<=k; i++)
        { 
			p=u*x1; 
			q=v*y1;
            pq=(u+v)*(x1+y1);
            u=p-q+ar[i]; 
			v=pq-p-q+ai[i];
        }
        
		g1=u*u+v*v;
        if (g1>=g)
        { 
			if (is!=0)
            { 
				it=1;
                g65c(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,&k,&is,&it);
                if (it==0) 
					goto l40;
            }
            else
            { 
				g60c(&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;
                    g65c(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,&k,&is,&it);
                    if (it==0) 
						goto l40;
                }
            }
            
			g90c(xr,xi,ar,ai,&x,&y,&p,&w,&k);
        }
        else
        { 
			g=g1; 
			x=x1; 
			y=y1; 
			is=0;
            if (g<=1.0e-22)
				g90c(xr,xi,ar,ai,&x,&y,&p,&w,&k);
            else
            { 
				u1=k*ar[0]; 
				v1=ai[0];
                for (i=2; i<=k; i++)
                { 
					p=u1*x; 
					q=v1*y; 
					pq=(u1+v1)*(x+y);
                    u1=p-q+(k-i+1)*ar[i-1];
                    v1=pq-p-q+(k-i+1)*ai[i-1];
                }
                
				p=u1*u1+v1*v1;
                if (p<=1.0e-20)
                { 
					it=0;
                    g65c(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,&k,&is,&it);
                    if (it==0) 
						goto l40;
                    
					g90c(xr,xi,ar,ai,&x,&y,&p,&w,&k);
                }
                else
                { 
					dx=(u*u1+v*v1)/p;
                    dy=(u1*v-v1*u)/p;
                    t=1.0+4.0/k;
                    g60c(&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;
                        g65c(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,&k,&is,&it);
                        if (it==0) 
							goto l40;
                    }
                    
					g90c(xr,xi,ar,ai,&x,&y,&p,&w,&k);
                }
            }
        }
        
		if (k==1) 
			jt=0;
        else 
			jt=1;
    }
    
	return TRUE;
}

//////////////////////////////////////////////////////////////////////
// 内部函数
//////////////////////////////////////////////////////////////////////
void CNLequation::g60c(double* t,double* x,double* y,double* x1,double* y1,double* dx,double* dy,double* p,double* q,int* k,int* it)
{ 
	*it=1;
    while (*it==1)
    { 
		*t=*t/1.67; 
		*it=0;
        *x1=*x-*t*(*dx);
        *y1=*y-*t*(*dy);
        if (*k>=30)
		{ 
			*p=sqrt(*x1*(*x1)+*y1*(*y1));
            *q=exp(75.0/(*k));
            if (*p>=*q) 
				*it=1;
        }
    }
}

//////////////////////////////////////////////////////////////////////
// 内部函数
//////////////////////////////////////////////////////////////////////
void CNLequation::g90c(double xr[],double xi[],double ar[],double ai[],double* x,double* y,double* p,double* w,int* k)
{ 
	int i;
    for (i=1; i<=*k; i++)
    { 
		ar[i]=ar[i]+ar[i-1]*(*x)-ai[i-1]*(*y);
        ai[i]=ai[i]+ar[i-1]*(*y)+ai[i-1]*(*x);
    }
    
	xr[*k-1]=*x*(*w); 
	xi[*k-1]=*y*(*w);
    *k=*k-1;
    if (*k==1)
    { 
		*p=ar[0]*ar[0]+ai[0]*ai[0];
        xr[0]=-*w*(ar[0]*ar[1]+ai[0]*ai[1])/(*p);
        xi[0]=*w*(ar[1]*ai[0]-ar[0]*ai[1])/(*p);
    }
}

//////////////////////////////////////////////////////////////////////
// 内部函数
//////////////////////////////////////////////////////////////////////
void CNLequation::g65c(double* x,double* y,double* x1,double* y1,double* dx,double* dy,double* dd,double* dc,double* c,int* k,int* is,int* it)
{ 
	if (*it==0)
    { 
		*is=1;
        *dd=sqrt(*dx*(*dx)+*dy*(*dy));
        if (*dd>1.0) 
			*dd=1.0;
        *dc=6.28/(4.5*(*k)); 
		*c=0.0;
    }
    
	while(TRUE)
    { 
		*c=*c+*dc;
        *dx=*dd*cos(*c); 
		*dy=*dd*sin(*c);
        *x1=*x+*dx; 
		*y1=*y+*dy;
        if (*c<=6.29)
        { 
			*it=0; 
			return;
		}
        
		*dd=*dd/1.67;
        if (*dd<=1.0e-07)
        { 
			*it=1; 
			return;
		}
        
		*c=0.0;
    }
}

//////////////////////////////////////////////////////////////////////
// 求非线性方程一个实根的蒙特卡洛法
//
// 调用时，须覆盖计算方程左端函数f(x)值的虚函数double Func(double x)
//
// 参数：
// 1. double* x - 传入初值（猜测解），返回求得的实根
// 2. double xStart - 均匀分布的端点初值
// 3. int nControlB - 控制参数
// 4. double eps - 控制精度，默认值为0.000001
//
// 返回值：无
//////////////////////////////////////////////////////////////////////
void CNLequation::GetRootMonteCarlo(double* x, double xStart, int nControlB, double eps /*= 0.000001*/)
{ 
	int k;
    double xx,a,y,x1,y1,r;
    
	// 求解条件
	a = xStart; 
	k = 1; 
	r = 1.0;

	// 初值
	xx = *x; 
	y = Func(xx);

	// 精度控制求解
    while (a>=eps)
    { 
		x1=rnd(&r);

		x1=-a+2.0*a*x1;
        x1=xx+x1; 
		y1=Func(x1);
        
		k=k+1;
        if (fabs(y1)>=fabs(y))
        { 
			if (k>nControlB) 
			{ 
				k=1; 
				a=a/2.0; 
			}
		}
        else
        { 
			k=1; 
			xx=x1; 
			y=y1;
            if (fabs(y)<eps)
            { 
				*x=xx; 
				return;