Newton Raphson 混合算法无法解决

Posted 2023-03-12

【中文标题】Newton Raphson 混合算法无法解决

【英文标题】：Newton Raphson hybrid algorithm not reaching a solution

【发布时间】：2012-11-02 00:50:52

【问题描述】：

问题的简要说明：我使用 Newton Raphson 算法在多项式中求根，但在某些情况下不起作用。为什么？

我从“c++ 中的数值配方”中提取了一个 Newton 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 个）其他非常相似的多项式（精度为 1e-10 且迭代次数很少！）时，这个多项式没有达到准确的错误！

我知道这个问题很难，而且可能没有真正直接的答案，但我被这个问题困扰了几天，我不知道如何解决它。非常感谢您花时间阅读我的问题。

【问题讨论】：

+1 问题很清楚，代码在那里，还有多余的漂亮图。

@AlessandroTeruzzi 图片在那里，因此可以看出多项式没有零导数或奇怪的行为，这是 New-raph 的常见问题，但如果它是多余的，我可以退出该图。代码可能会出现在那里和那里，但它取自剑桥大学制作的 C++ 书中的数字食谱。它与书中发现的几乎没有任何变化，这本书在谈论数值方法时非常值得信赖。顺便说一句，感谢您的 +1，这对吸引可能知道如何解决我的问题的人很有帮助。

您能否展示您的最小 C++ 程序，该程序使用给定的多项式调用此函数。我猜应该是 10-15 行代码。

@JonathanLeffler 我可以这样做，如果它会澄清更多，但除了我已经写过的内容之外，您唯一必须知道的是我将数组中多项式的系数从低到高排序。我构造了一个 double* coeffs=new double[degree];数组，从低度到高度系数填充它并调用函数。其他变量是度数（在本例中为 4）和我已经发布的 4 个变量。如果您仍然觉得有必要举例说明，我会毫无问题地发布一个。

这看起来像是一个愚蠢的建议，但您是否考虑过在代码中添加日志，可能会帮助您查看问题所在...

【参考方案1】：

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

   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为负数。很有可能你是对的；浮点运算的相等性很少是一个好主意。

嗨，@JonathanLeffler 和 sizzzzlerz。实际上，这不是代码中返回的常见情况。仅当结果恰好在括号的角上时，该返回才有效。如果您检查下面的某些行，则写了一个“if”表达式，这实际上就是您向我提出的建议。那是代码在 99% 的时间结束的返回行。如前所述，我已经尝试过包含 1000 多个类似多项式的代码，并且它适用于 1e-10 的精度。

再仔细看代码，发现if (xl == rts)条件是测试一个加法后候选根是否发生了变化；我认为这是一个合理的比较（尽管必须小心；编译器可能会把事情搞砸，但可能不会）。主要收敛标准是if (fabs(dx) < xacc) return rts;// Convergence criterion，这是合理制定的（它查看 Δx 是否小于请求的限制）。

@JonathanLeffler 是的！但在这种情况下，fabs(dx)=0.15，在 1000 次迭代后我发现这个值非常高

【参考方案2】：

我不确定具体原因，但在这种情况下，检查函数是否足够快地减小似乎不起作用。

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

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

【讨论】：

嗯，很有趣....明天我会检查它并告诉你答案。谢谢

@AnderBiguri：我想我在你的代码中发现了一个问题——我更新了我的答案。

没错！哇，我已经花了这么多小时来审查代码，却没有注意到那个小（同时也是巨大的！）错误。非常感谢！