Newton Raphson hybrid algorithm not reaching a solution

Question

Brief explanation of the problem: I use Newton Raphson algorithm for root finding in polynomials and doesn't work in some cases. why?

I took from "numerical recipes in c++" a Newton Raphson hybrid algorithm, which bisects in case New-Raph is not converging properly (with a low derivative value or if the convergence speed is not fast).

I checked the algorithm with several polynomials and it worked. Now I am testing in inside the software I have and I always got an error with an specific polynomial. My problem is that I don't know why this polynomial just doesn't get to the result, when much others do. As I want to improve the algorithm for any polynomial y need to know which one is the reason of no convergence so I can treat it properly.

Following I will post all the information I can provide about the algorithm and the polynomial in which I have the error.

The polynomial:

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

It's first derivative:

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

Plot: enter image description here

Roots (by Matlab):

  -2.133112008595826          1.371976341295347          0.883715461977390 
  -0.679837109933505

Algorithm:

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.
}

The algorithm is called with the next variables:

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;

the problem is that the algorithm exceeds the amximum iterations and there is an arror to the root of aproximatedly 0.15.

So my direct and abrebiated question is: Why this polynomial does not reach an accurate error when many (1000 at least) other very similar polynomials do (wuth 1e-10 of precision and few iterations!)

I know the question is difficult and it may not have a really direct answer, but I am stuck with this for some days and I don't know how to solve it. Thank you very much for taking time for reading my question.

+1 the question is clear, the code is there and there is also a redundant pretty picture. — Alessandro Teruzzi, Nov 13 '12 at 15:34
@AlessandroTeruzzi The picture is there so it can be seen that the polynomial has no zero derivatives or strange behaviour, a usual problem with New-raph, but I can quit the figure if it is redundant. The code may appear there and there, but it is taken from Numerical Recipes in c++ book, made in Cambridge university. It has almost none change from the one found in the book, a book pretty trustable when talking about numerical methods. Thanks for the +1 by the way, it helps a lot for taking attention for the people that may know how to solve my problem. — Ander Biguri, Nov 13 '12 at 15:45
Can you show your minimal C++ program that invokes this function with the given polynomial. I'm guessing that should be 10-15 lines of code. — Jonathan Leffler, Nov 13 '12 at 16:01
@JonathanLeffler I can do that if it will clarify more, but the only thing you must know besides what I already wrote is that I order the coefficients of the polynomials in the array from low to high degree. I construct a double* coeffs=new double[degree]; array, fill it from low degree to high degree coefficients and call the function. The other variables are the degree (in this case 4) and the 4 variables I already posted. If you still find the necessity of an example say it and I will post one without any problem. — Ander Biguri, Nov 13 '12 at 16:06
this will look like a stupid suggestion, but did you think about adding a log to your code, might help you to see what is the problem... — NoSenseEtAl, Nov 13 '12 at 16:33
@NoSenseEtAl sorry, I am half-new to programming, what do you mean with a log? Go saving what the variables are in each iteration? — Ander Biguri, Nov 13 '12 at 16:48
yes, you can like cout values for the first 100 iterations, "real logging" has has levels of logging (INFO, WARNING... ) if you want to know more in general check out google glog lib. — NoSenseEtAl, Nov 13 '12 at 17:30

score 3 · Answer 1 · answered Nov 13 '12 at 15:28

3

Without having run your code, my initial guess is that you are comparing to floating point values for equality to determine if your solution has converged.

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

Maybe you should calculate it as a ratio:

   diff = fabs(xl - rts);
   if (diff/xl <= 1.0e-8)  // pick your own accuracy value here
      return rts;

answered Nov 13 '12 at 15:28

sizzzzlerz

3,905
3
20
33

Probably need `fabs(xl)` in the divisor, in case `xl` is negative. There's a decent chance you're right; equality for floating point arithmetic is seldom a good idea. – Jonathan Leffler Nov 13 '12 at 15:49
Hi,@JonathanLeffler and sizzzzlerz. Actually that is not the usual case for the return in the code. That return will only work if the result is exactly in the corner of the brackets. If you check some lines below there written an 'if' expression that is practically what you propose me. That one is the return line where the code ends the 99% of the times. As said before, I have tried the code with over 1000 similar polynomials and it works with 1e-10 of precission. – Ander Biguri Nov 13 '12 at 15:52
1

Looking at the code harder, I see that the `if (xl == rts)` condition is testing whether the candidate root has changed after an addition; that is a legitimate comparison, I think (though one has to be careful; and the compilers might mess things up, but probably won't). The main convergence criterion is `if (fabs(dx) < xacc) return rts;// Convergence criterion`, and that is reasonably formulated (it looks to see if the Δx is smaller than the limit requested). – Jonathan Leffler Nov 13 '12 at 16:00
@JonathanLeffler yes! but in this case fabs(dx)=0.15, something that I found really high after 1000 iterations – Ander Biguri Nov 13 '12 at 16:10

Vaughn Cato · Accepted Answer · 2012-11-13T18:13:16.523

3

I'm not sure exactly why, but the check to see if the function is decreasing fast enough doesn't appear to work in this case.

It works if I do it like this:

  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;

UPDATE

It looks like there is a problem in your code:

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

should be

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

edited Nov 13 '12 at 18:13

answered Nov 13 '12 at 16:23

Vaughn Cato

59,967
5
75
116

Hummm interesting.... I will check it tomorrow and let you know the answer. Thanks – Ander Biguri Nov 13 '12 at 16:25
@AnderBiguri: I think I found a problem in your code -- I updated my answer. – Vaughn Cato Nov 13 '12 at 18:13
Exactly! Wow I've been so many hours reviewing the code without noticing that small (and huge in the same time!) error. thank you very much! – Ander Biguri Nov 14 '12 at 10:07

Newton Raphson hybrid algorithm not reaching a solution

2 Answers2