nonlinearequation.inl

vc++实现线性方程组求解 1全选主元高斯消元法 2全选主元高斯-约当消元法 3三对角方程组的追赶法 4一般带型方程组求解 5对称方程组的分解法 6对称正定方程组的平方根法

//NonLinearEquation.inl 	非线性方程(组)求解函数(方法)定义
// Ver 1.0.0.0
// 版权所有(C) 何渝, 2002
// 最后修改: 2002.5.31.

#ifndef _NONLINEAREQUATION_INL
#define _NONLINEAREQUATION_INL

//用二分法搜索方程f(x)=0在区间[a,b]内的全部实根
template <class _Ty>
inline size_t 
RootHalves(_Ty a, _Ty b, _Ty step, _Ty eps, valarray<_Ty>& x, size_t m)
{
    int n(0), js;
    _Ty z=a, z1, y, y1, z0, y0;
   
	y = FunctionValueRH(z);
    while((z<=b+step/2.0)&&(n!=m))
	{ 
	  if(Abs(y) < eps)
      {
		n = n + 1;
		x[n-1] = z;
		z = z + step / 2.0;
		y=FunctionValueRH(z);
      }
      else
      {
		z1=z+step;
		y1=FunctionValueRH(z1);
        if(Abs(y1)<eps)
        {
			n=n+1;
			x[n-1]=z1;
            z=z1+step/2.0;
			y=FunctionValueRH(z);
        }
        else
			if(y*y1>0.0)
            {
				y=y1;
				z=z1;
			}
            else
            {
				js=0;
                while(js==0)
                {
					if(Abs(z1-z)<eps)
                    {
						n=n+1;
						x[n-1]=(z1+z)/2.0;
                        z=z1+step/2.0;
						y=FunctionValueRH(z);
                        js=1;
                    }
                    else
                    {
						z0=(z1+z)/2.0;
						y0=FunctionValueRH(z0);
                        if(Abs(y0)<eps)
                        {
							x[n]=z0; n=n+1; js=1;
                            z=z0+step/2.0;
							y=FunctionValueRH(z);
                        }
                        else
							if((y*y0) < 0.0)
							{
								z1 = z0;
								y1 = y0;
							}
							else
							{
								z = z0;
								y = y0;
							}
                      }
                  }
              }
          }
	}
    return(n);	//返回搜索到的实根个数
}

//牛顿(Newton)法求解非线性方程一个实根
template <class _Ty>
inline int 
RootNewton( _Ty& x, _Ty eps, size_t js)
{ 
    int k = 1, r = js;
    _Ty  x1;
	_Ty x0 = x;
	valarray<_Ty> y(2);

    FunctionValueRN(x0, y);
    _Ty d=eps + 1.0;
    while((d > eps || FloatEqual(d, eps))&&(r != 0))
    {
		if(FloatEqual(y[1],0))	return(-1);
        x1 = x0 - y[0] / y[1];
        FunctionValueRN(x1,y);
        d=Abs(x1 - x0);
		_Ty p=Abs(y[0]);
        if(p > d) d = p;
        x0 = x1;
		r -= 1;
    }
    x = x1;
    k = js - r;
    return(k);
}

//埃特金(Aitken)法求解非线性方程一个实根
template <class _Ty>
inline int 
RootAitken( _Ty& x, _Ty eps, size_t js)
{
    _Ty u, v, x0;
    int r(0), flag(0);

	x0 = x;

    while((flag == 0) && (r != js))
    {	
		r++; 
        u = FunctionValueRA(x0);
		v = FunctionValueRA(u);
        if(Abs(u-v)<eps)
		{
			x0 = v;
			flag = 1;
		}
        else x0=v-(v-u)*(v-u)/(v-2.0*u+x0);
    }
    x = x0;
	r = js - r;
    return(r);
}

//连分式(Fraction)法求解非线性方程一个实根
template <class _Ty>
inline int 
RootFraction(_Ty& x, _Ty eps)
{
	int r(10), i, j, m, it;
    _Ty x0, q(1.0e+35), h(0), z;
	valarray<_Ty> a(10), y(10);
	
	x0 = x; 

    while(r!=0)
    { 
		r = r - 1; 
		j = 0; 
		it = r;
        while(j<=7)
        {
			if(j<=2)
				z = x0 + 0.1 * j;
            else
				z = h;
            y[j] = FunctionValueRF(z);
            h = z;
            if(j==0)
				a[0] = z;
            else
            {
				m = 0;
				i = 0;
                while((m==0)&&(i<j))
                {
					if(FloatEqual((h-a[i]),0))
						m = 1;
                    else
						h=(y[j]-y[i])/(h-a[i]);
                    i++;
                }
                a[j] = h;
                if(m!=0) a[j]=q;
                h=0.0;
                for(i=j-1; i>=0; i--)
                {
					if(FloatEqual((a[i+1]+h),0))
						h = q;
                    else
						h=-y[i]/(a[i+1]+h);
                }
                h = h + a[0];
            }
            if(Abs(y[j])>=eps) j++;
            else 
			{
				j = 10;
				r = 0;
			}
        }
        x0 = h;
    }
    x = h;
    return(10-it);
}

//QR法求代数方程全部根
template <class _Ty, class _Tz>
inline int 
RootQR(valarray<_Ty>& a, valarray<_Tz>& x, _Ty eps, size_t jt)
{
	int i, j, n;

    n = x.size();			//多项式方程的次数

	matrix<_Ty> q(n, n);

    for(j=0; j<n; j++)
		q(0,j) = -a[n-j-1] / a[n];

    for(j=1; j<n; j++)
		for(i=0; i<n; i++)
			q(j,i) = 0;

	for(i=0; i<n-1; i++)
      q(i+1, i) = 1;
	
	i = EigenvalueVectorHessenbergQR(q,x,eps,jt);	//求全部特征值的QR方法
    
	return(i);
}

//牛顿下山(NewtonHillDown)法求解实系数代数方程全部根(实根和复根)
template <class _Ty, class _Tz>
inline int 
RootNewtonHillDown(valarray<_Ty>& a, valarray<_Tz>& cx)
{
	int m, i, jt, k, is, it, n;
    _Ty t, x, y, x1, y1, dx, dy, dd, dc, c, g, g1;
	_Ty w, u, p, q, pq, v, u1, v1;
  
    n = cx.size();			//多项式方程的次数
    m = n;

    while((m > 0) && (FloatEqual(a[m], 0.0))) m--;

    if(m <= 0) return(-1);
    
	for(i = 0; i <= m; i++) a[i] /= a[m];

    for(i = 0; i <= m / 2; i++)
    {
		w = a[i];
		a[i] = a[m - i];
		a[m - i] = w;
	}

    k = m;
	is = 0;
	w = 1.0;
    jt = 1;
    while(jt == 1)
    {
		pq = Abs(a[k]);
		while(pq < LONGDOUBLEERROR)
        {
			cx[k - 1] = complex<_Ty>(0.0, 0.0);
			k = k - 1;
            if(k == 1)
            {
				cx[0] = complex<_Ty>(-a[1] * w / a[0], 0.0);
                return(1);
            }
            pq = Abs(a[k]);
        }
		q = log(pq);
		q = q / (1.0 * k);
		q = exp(q);
        p = q;
		w = w * p;
        for(i = 1; i <= k; i++)
        {
			a[i] /= q;
			q *= p;
		}
        x = 0.0001;
		x1 = x; 
		y = 0.2;
		y1 = y; 
		dx = 1.0;
        g = 1.0e+37; 
 
zjn:    u = a[0];
		v = 0.0;
        for(i = 1; i <= k; i++)
        { 
			p = u * x1;
			q = v * y1;
            pq = (u + v) * (x1 + y1);
            u = p - q + a[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 zjn;
			}
			else
			{
				g60(t, x, y, x1, y1, dx, dy, p, q, k, it);
				if(t >= 1.0e-03)  goto zjn;
				if(g > 1.0e-18)
				{
					it = 0;
					g65(x, y, x1, y1, dx, dy, dd, dc, c, k, is, it);
					if(it == 0) goto zjn;
				}

			}
			g90(cx , a, x, y, p, q, w, k);
        }
        else
        {
			g = g1;
			x = x1;
			y = y1;
			is = 0;
            if(g <= 1.0e-22)	
				g90(cx, a, x, y, p, q, w, k);
            else
            {
				u1 = k * a[0]; 
				v1 = 0.0;
                for(i = 2; i <= k; i++)
                {
					p = u1 * x;
					q = v1 * y;
					pq = (u1 + v1) * (x + y);
                    u1 = p - q + (k - i + 1) * a[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 zjn;
                    g90(cx, a, 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 zjn;
                    if(g > 1.0e-18)
                    {
						it = 0;
                        g65(x, y, x1, y1, dx, dy, dd, dc, c, k, is, it);
                        if(it == 0) goto zjn;
                    }
                    g90(cx, a, x, y, p, q, w, k);
                }
            }
        }

        if(k == 1) jt = 0;
        else jt = 1;
    }
    return(1);
}

//牛顿下山(NewtonHillDown)法求解复系数代数方程全部根(实根和复根)
//重载RootNewtonHillDown()
template <class _Ty>
inline int 
RootNewtonHillDown(complex<_Ty> a[], valarray<complex<_Ty> >& cx)
{
	int m, i, jt, k, is, it, n;
    _Ty t, x, y, x1, y1, dx, dy, dd, dc, c, g, g1, mo, sb, xb;
	_Ty w, u, p, q, pq, v, u1, v1;
    
	n = cx.size();			//多项式方程的次数
    m = n;
	
    while((m > 0) && (FloatEqual(Abs(a[m]), 0))) m--;

    if(m <= 0) return(-1);

    mo = Abs(a[m]);

	for(i = 0; i <= m; i++) a[i] /= mo;

    for(i = 0; i <= m / 2; i++)
    {
		swap(a[i], a[m-i]);		
	}

    k = m;
	is = 0;
	w = 1.0;
    jt = 1;
    while(jt == 1)
    {
		pq = Abs(a[k]);
		while(pq < LONGDOUBLEERROR)
        {
			cx[k - 1] = complex<_Ty>(0.0, 0.0);
			k--;
            if(k == 1)
            {
				mo = a[0].real() * a[0].real() + a[0].imag() * a[0].imag();
				sb = -w * (a[0].real() * a[1].real() + a[0].imag() * a[1].imag()) / mo;
				xb =  w * (a[1].real() * a[0].real() - a[0].imag() * a[1].imag()) / mo;
				cx[0] = complex<_Ty>(sb, xb);
                return(1); 
            }
			pq = Abs(a[k]);
        }

		q = log(pq);
		q = q / (1.0 * k);
		q = exp(q);
        p = q;
		w = w * p;
        
		for(i = 1; i <= k; i++)
        {
			a[i] /= q;
			q *= p;
		}
        
		x = 0.0001;
		x1 = x; 
		y = 0.2;
		y1 = y; 
		dx = 1.0;
        g = 1.0e+37; 
 
zjn:   
		u = a[0].real();
		v = a[0].imag();
        
		for(i = 1; i <= k; i ++)
        { 
			p = u * x1;
			q = v * y1;
            pq = (u + v) * (x1 + y1);
            u = p - q + a[i].real();
			v = pq - p - q + a[i].imag();
        }
        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 zjn;
			}
			else
			{
				gg60(t, x, y, x1, y1, dx, dy, p, q, k, it);
				if(t >= 1.0e-03)  goto zjn;
				if(g > 1.0e-18)
				{
					it = 0;
					g65(x, y, x1, y1, dx, dy, dd, dc, c, k, is, it);
					if(it == 0) goto zjn;
				}
			}

			gg90(cx, a, x, y, p, q, w, k);
        }
        else
        {
			g = g1;
			x = x1;
			y = y1;
			is = 0;

            if(g <= 1.0e-22)	
				gg90(cx, a, x, y, p, q, w, k);
            else
            {
				u1 = k * a[0].real(); 
				v1 = a[0].imag();
                for(i = 2; i <= k; i++)
                {
					p = u1 * x;
					q = v1 * y;
					pq = (u1 + v1) * (x + y);
                    u1 = p - q + (k - i + 1) * a[i - 1].real();
                    v1 = pq - p - q + (k - i + 1) * a[i - 1].imag();
                }
                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 zjn;
                    gg90(cx, a, 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;
                    gg60(t, x, y, x1, y1, dx, dy, p, q, k, it);
                    if(t >= 1.0e-03) goto zjn;
                    if(g > 1.0e-18)