Search code examples
cwhile-loopnewtons-method

Newton convergence method not working


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.


Solution

  • 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>