牛顿-拉夫逊混合算法无法达到解决方案

9
问题简要说明:我在多项式求根中使用了牛顿-拉夫森算法,在某些情况下无法正常工作。为什么?
我从《C++数值计算》中借鉴了一个牛顿-拉夫森混合算法,如果牛顿-拉夫森算法不收敛(导数值低或收敛速度不快),则采用二分法。
我检查了该算法在多个多项式上的表现,并且它可以正常工作。现在我正在软件中测试它,在特定的多项式中总是出现错误。我的问题是,我不知道为什么这个多项式不能得到结果,而其他多项式可以。因为我想改进该算法以适用于任何多项式,所以我需要知道哪个多项式导致不收敛,以便我可以正确地处理它。
接下来,我将发布有关该算法和多项式错误的所有信息。
多项式:
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个)非常相似的多项式达到了1e-10的精度和少量迭代,但是这个多项式无法达到精确的误差?
我知道这个问题很难,可能没有一个真正的直接答案,但我被困扰了几天,不知道该怎么解决。非常感谢您花时间阅读我的问题。

3
问题清晰,代码已经存在,还有一个多余的漂亮图片。+1 - Alessandro Teruzzi
你能展示一下调用这个给定多项式函数的最小C++程序吗?我猜应该是10-15行代码。 - Jonathan Leffler
@JonathanLeffler 如果这样做可以更清晰明了,那我可以这样做。但是除了我已经写的内容之外,你必须知道的唯一一件事情就是我按照多项式系数的次数从低到高的顺序将它们存储在数组中。我构造了一个double* coeffs=new double[degree]; 数组,从低次到高次系数填充它,并调用函数。其他变量是度数(在这种情况下为4)和我已经发布的4个变量。如果你仍然觉得需要一个示例,请告诉我,我会毫不费力地发布一个。 - Ander Biguri
2
这可能听起来很愚蠢,但你有没有考虑在你的代码中添加日志,这可能会帮助你看到问题所在... - NoSenseEtAl
是的,您可以像 cout 值一样输出前 100 次迭代,“真实日志”具有日志级别(INFO、WARNING...),如果您想了解更多,请查看 Google glog 库。 - NoSenseEtAl
显示剩余2条评论
2个回答

3

在没有运行您的代码之前,我最初的猜测是,您正在比较浮点值的相等性以确定您的解决方案是否已经收敛。

   if (xl == rts) 
        return rts; //Change in root is negligible.

也许你应该将它计算为比率:
   diff = fabs(xl - rts);
   if (diff/xl <= 1.0e-8)  // pick your own accuracy value here
      return rts;

可能需要在除数中加入 fabs(xl),以防 xl 为负数。你说的很有可能是对的;对于浮点运算来说,相等很少是一个好主意。 - Jonathan Leffler
嗨,@JonathanLeffler和sizzzzlerz。实际上,这不是代码中通常的返回情况。如果结果恰好在括号的角落里,那么该返回语句才能起作用。如果您向下检查几行,会发现一个“if”表达式,这实际上就是您建议我的内容。这是代码99%的结束行。如前所述,我已经尝试了1000多个类似的多项式代码,并且它可以精确到1e-10。 - Ander Biguri
1
仔细查看代码,我发现 if (xl == rts) 条件是在测试候选根是否在加法后发生了变化;我认为这是一个合理的比较(尽管需要小心;编译器可能会搞砸,但可能不会)。主要的收敛准则是 if (fabs(dx) < xacc) return rts;// 收敛准则,这是一个合理的公式(它检查 Δx 是否小于所请求的限制)。 - Jonathan Leffler
@JonathanLeffler 是的!但在这种情况下,fabs(dx)=0.15,在1000次迭代后我发现这个值非常高。 - Ander Biguri

3

我不确定为什么,但检查函数是否下降得够快似乎在这种情况下无效。

如果像下面这样操作,它是有效的:

  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);

嗯,有趣……我明天会检查一下并告诉你答案。谢谢。 - Ander Biguri
@AnderBiguri:我认为我在你的代码中发现了一个问题——我更新了我的答案。 - Vaughn Cato
没错!哇,我花了这么多时间来审查代码,竟然没有注意到那个小错误(同时又是一个巨大的错误!)。非常感谢! - Ander Biguri

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接