Newton convergence method not working

Question

I am trying to approximate the root of a polynomial using Newton-Raphson method. The code I wrote to do it looks like this:

#include <stdio.h>
#include <math.h>

int main (){

double c, nq, nnq, nnnq, v, h, q, o;

o=2;
c=-0.55;
nq=-0.04625;
nnq=-0.55;
nnnq=1;

   while(fabs(c-v) > 0.000001)
   {
      nq=((2*(o-1)+1)*(c)*(nnq)-(o-1)*nnnq)/((o-1)+1);                                                     
      q=(-o*c*nq+o*nnq)/(1-(c*c));
      h=(c-(nq/q));
      printf("The better root is %.15lf\n",h);
      v=c;
      c=h;
   }

}

I know it is not necessary to write the variables o,c,nq, etc sin I could just use their exact values. This is a part of a larger problem and I need those variables, so ignore that.

This program output this:

The better root is -0.578030303030303
The better root is -0.591696792857493
The better root is -0.598609887802599
The better root is -0.602171714355970
The better root is -0.604024260228500
The better root is -0.604992519745332
The better root is -0.605499890229896
The better root is -0.605766110042157
The better root is -0.605905895095070
The better root is -0.605979319651017
The better root is -0.606017894664121
The better root is -0.606038162857992
The better root is -0.606048812800124
The better root is -0.606054408979837
The better root is -0.606057349623975
The better root is -0.606058894866533
The better root is -0.606059706860161

When instead it should converge to the point -0.57735026918963. I know the Newton-Raphson converges for sure, so the error should be on the code. I've also tried to localitzate the problem using printf, and I think The problem comes in the second iteration. I think the program fails to calculate nq correctly but I don't know why.

One serious problem (although not the major one) is that `v` is uninitialised. This leads me to believe that you are compiling without warnings enabled (or perhaps just ignoring warnings) ? — Paul R, May 14 '18 at 18:25
can you add the polynomial that you want to calculate to the question too? — Afshin, May 14 '18 at 18:29
@PaulR True, in the other program I have v=-20 to make sure the while loop starts. Either way, the output is the same. — John Keeper, May 14 '18 at 18:33
@Afshin It is the [Legendre polynomial](https://en.wikipedia.org/wiki/Legendre_polynomials) for n=2: (1/2)*(3*x^2-1) — John Keeper, May 14 '18 at 18:34
Are you sure that first derivative term is correct? Looks funny to me. — Bathsheba, May 14 '18 at 18:35
@Bathsheba Yes, to calculate the derivatives of Legendre polynomials, you can use [this](https://i.gyazo.com/b2d8640af132bbf9e50dd94a411e91e4.png) recursive method, knowing that p_0=1 and p_1=x. In my case, n=2. — John Keeper, May 14 '18 at 18:37
@Bathsheba yea,it looked funny to me too. That's why I asked for equation. — Afshin, May 14 '18 at 18:37
@JohnKeeper then you are using 2 different recursive equation at same time. I think that's why you cannot converge. Try to use newton method only on your equation (manually add differential to code). I guess you should converge that time. — Afshin, May 14 '18 at 18:39
Why did you write `((o-1)+1)` instead of `o`? Perhaps the parenthesis is misplaced. — rici, May 14 '18 at 18:49

Afshin · Answer 1 · 2018-05-14T19:04:00.800

1

This is newton method for your equation(it is a quick code, don't check variable name):

#include <stdio.h>
#include <math.h>

int main ()
{    
    double s = 2.0, fx = 0, dfx = 0, p = 0;

    while(fabs(s - p) > 0.000001)
    {
        fx = 0.5 * (3 * s * s - 1);
        dfx = 3 * s;
        p = s;
        s = s - (fx / dfx);
        printf("The better root is %.15lf\n", s);
    }    
    return 0;
}

and it converged to 0.577350269189626. Your problem is that you are trying to calculate 2 recursions at same time. Btw, in your question you said that you want to calculate "root of a polynomial". I didn't get exactly what you meant. If from root you mean square root of your equation, you need to update this code and change fx and dfx accordingly.

edited May 14 '18 at 19:04

answered May 14 '18 at 18:58

Afshin

8,839
1
18
53

`> 0.000001` in reality would need more thought. – Bathsheba May 14 '18 at 19:38
@Bathsheba yea,I just copy/paste his code and edited it. :) in addition start point is important too (I started at 2.0). – Afshin May 14 '18 at 19:42

Lutz Lehmann · Accepted Answer · 2018-05-22T18:10:41.110

Instead of just computing x = sqrt(1.0/3), you want to incorporate the recursive evaluation of the Legendre polynomials up to order o, probably to extend the method later to values of o larger than 2. The iteration is

P(0,c) = 1; P(1,c) = c;
(n+1)*P(n+1,c) = (2*n+1)*c*P(n,c) - n*P(n-1,c),   n=1, ... , o-1

and the derivative can be computed as

(1-c^2)*P'(o,c) = -n*c*P(o,x) + n*P(o-1,c).

This iteration you would need to include inside the Newton loop, in the ideal case using an object for the Legendre polynomial with methods for value and derivative. I've modified your structure to work in JavaScript:

var my_div = document.getElementById("my_div");
var c = -0.55;
var v = -20;
var o = 2;
while( Math.abs(v-c) > 1e-12 ) {
    p0 = 1; 
    p1 = c;
    for(n=1; n<o; n++) {
        // p0, p1, p2 stand for P(n-1,c), P(n,c), P(n+1,c)
        p2 = ((2*n+1)*c*p1 - n*p0) / (n+1)
        p0 = p1; p1 = p2;
    }
    // p0, p1 now contain p(o-1,x), p(o,x)
    p1prime = ( -o*c*p1 + o*p0) / (1-c*c);
    h = c - p1/p1prime;
    my_div.innerHTML += "<p>The better root is "+h+"</p>";
    v = c; c = h;
 }

<div id="my_div"></div>

Newton convergence method not working

2 Answers2