Newton Raphson混合算法没有达到解决方案

Question

问题的简要说明:我使用Newton Raphson算法在多项式中进行根发现,但在某些情况下不起作用.为什么？

我从"c ++中的数字配方"中得到了一个N​​ewton Raphson混合算法,如果New-Raph没有正确收敛(具有低导数值或者如果收敛速度不快),则将其一分为二.

我用几个多项式检查了算法并且它有效.现在我在我的软件内部进行测试,并且我总是遇到一个特定多项式的错误.我的问题是,我不知道为什么这个多项式只是没有得到结果,而其他人做的很多.因为我想改进任何多项式的算法,需要知道哪一个是没有收敛的原因,所以我可以正确对待它.

下面我将发布我可以提供的有关算法和我遇到错误的多项式的所有信息.

多项式:

f(t)= t^4 + 0,557257315256597*t^3 - 3,68254086033178*t^2 +
+ 0,139389107255627*t + 1,75823776590795

它的第一个衍生物:

 f'(t)= 4*t^3 + 1.671771945769790*t^2 - 7.365081720663563*t + 0.139389107255627

情节:

根(由Matlab提供):

  -2.133112008595826          1.371976341295347          0.883715461977390 
  -0.679837109933505

算法:

double rtsafe(double* coeffs, int degree, double x1, double x2,double xacc,double xacc2)
    {
    int j;
    double df,dx,dxold,f,fh,fl;
    double temp,xh,xl,rts;
    double* dcoeffs=dvector(0,degree);
    for(int i=0;i<=degree;i++)
        dcoeffs[i]=0.0;
    PolyDeriv(coeffs,dcoeffs,degree);
    evalPoly(x1,coeffs,degree,&fl);
    evalPoly(x2,coeffs,degree,&fh);
    evalPoly(x2,dcoeffs,degree-1,&df);
if ((fl > 0.0 && fh > 0.0) || (fl < 0.0 && fh < 0.0))
    nrerror("Root must be bracketed in rtsafe");

if (fl == 0.0) return x1;
if (fh == 0.0) return x2;

if (fl < 0.0) { // Orient the search so that f(xl) < 0.
    xl=x1;
    xh=x2;
} else {
    xh=x1;
    xl=x2;
}
rts=0.5*(x1+x2);    //Initialize the guess for root,
dxold=fabs(x2-x1);  //the "stepsize before last,"
dx=dxold;           //and the last step

evalPoly(rts,coeffs,degree,&f);
evalPoly(rts,dcoeffs,degree-1,&dx);

for (j=1;j<=MAXIT;j++) { //Loop over allowed iterations

    if ((((rts-xh)*df-f)*((rts-xl)*df-f) > 0.0) //Bisect if Newton out of range,
            || (fabs(2.0*f) > fabs(dxold*df))) { //or not decreasing fast enough.
        dxold=dx;
        dx=0.5*(xh-xl);
        rts=xl+dx;
        if (xl == rts) 
            return rts; //Change in root is negligible.
    } else {// Newton step acceptable. Take it.
        dxold=dx;
        dx=f/df;
        temp=rts;
        rts -= dx;
        if (temp == rts)
            return rts;
    }
    if (fabs(dx) < xacc) 
        return rts;// Convergence criterion
    evalPoly(rts,coeffs,degree,&f);
    evalPoly(rts,dcoeffs,degree-1,&dx);
    //The one new function evaluation per iteration.
    if (f < 0.0) //Maintain the bracket on the root.
        xl=rts;
    else
        xh=rts;

}
//As the Accuracy asked to the algorithm is really high (but usually easily reached)
//the results precission is checked again, but with a less exigent result
dx=f/df;
if(fabs(dx)<xacc2)
    return rts;
nrerror("Maximum number of iterations exceeded in rtsafe");
return 0.0;// Never get here.
}

使用下一个变量调用该算法:

x1=0.019
x2=1.05
xacc=1e-10
xacc2=0.1
degree=4
MAXIT=1000
coeffs[0]=1.75823776590795;
coeffs[1]=0.139389107255627;
coeffs[2]=-3.68254086033178;
coeffs[3]=0.557257315256597;
coeffs[4]=1.0;

问题是该算法超过了最大迭代次数,并且以近似的方式存在一个错误0.15.

因此,我的直接和相关问题是:为什么当多个(至少1000个)其他非常相似的多项式执行时,这个多项式没有达到准确的误差(we 1e-10精度和几次迭代!)

我知道问题很难,而且可能没有真正直接的答案,但我坚持了几天,我不知道如何解决它.非常感谢你花时间阅读我的问题.

Answer 1

我不确定到底为什么，但检查函数是否下降得足够快在这种情况下似乎不起作用。

如果我这样做的话它会起作用：

  double old_f = f;

  .
  .
  .

    if ((((rts-xh)*df-f)*((rts-xl)*df-f) > 0.0) //Bisect if Newton out of range,
        || (fabs(2.0*f) > old_f)) { //or not decreasing fast enough.
  .
  .
  .

    if (fabs(dx) < xacc)
      return rts;// Convergence criterion
    old_f = f;

更新

您的代码似乎有问题：

evalPoly(rts,dcoeffs,degree-1,&dx);

应该

evalPoly(rts,dcoeffs,degree-1,&df);