rtsafe.cpp

Visual C++ 常用数值算法集 源代码

double rtsafe(double x1, double x2,double xacc)
{
	int maxit,j;
	double temp1,xl,xh,fl,fh,df,swap,temp,dxold,dum,f,dx;
	maxit = 100;
    funcd(x1, fl, df);
    funcd(x2, fh, df);
    if( (fl * fh) >= 0) cout<< "root must be bracketed"<<endl;
    if (fl < 0)
	{
        xl = x1;
        xh = x2;
	}
    else
	{
        xh = x1;
        xl = x2;
        swap = fl;
        fl = fh;
        fh = swap;
	}
    temp1 = 0.5 * (x1 + x2);
    dxold = fabs(x2 - x1);
    dx = dxold;
    funcd(temp1, f, df);
    for (j = 1; j<=maxit; j++)
	{
        dum = (temp1 - xh) * df - f * ((temp1- xl) * df - f);
        if ((dum >= 0)||( fabs(2 * f) > fabs(dxold * df)))
		{
            dxold = dx;
            dx = 0.5 * (xh - xl);
            temp1 = xl + dx;
            if (xl == temp1)  	
			{
				return temp1;
			    _c_exit();
			}
		}
        else
		{
            dxold = dx;
            dx = f / df;
            temp = temp1;
            temp1 = temp1 - dx;
            if (temp == temp1) 	
			{
				return temp1;
				_c_exit();
			}
        }
       if (fabs(dx) < xacc) 	
	   {
		   return temp1;
	      _c_exit();
	   }
       funcd(temp1, f, df);
        if (f < 0 )
		{
            xl = temp1;
            fl = f;
		}
        else
		{
            xh = temp1;
            fh = f;
        }
    }
	return temp1;
    cout<< "rtsafe exceeding maximum iterations."<<endl;
}