1

I've tried to implement Newton's method for polynomials. Like:

double xn=x0;
double gxn=g(w, n, xn);

int i=0;
while(abs(gxn)>e && i<100){
    xn=xn-(gxn/dg(w, n, xn));
    gxn=g(w, n, xn);

    i++;
}

where g(w, n, xn) computes the value of the function and dg(w, n, xn) computes the derivative.

As x0 I use starting point M which I found using Sturm's theorem.

My problem is that this method is divergent for some polynomials like x^4+2x^3+2x^2+2x+1. Maybe it's not regular, but I noticed that it happens when the solution of the equation is a negative number. Where can I look for an explanation?

Edit: dg

double result=0;
for(int i=0; i<n+1; i++)
    result+=w[i]*(n-i)*pow(x, n-i-1);

where n is the degree of polynomial

adolzi
  • 671
  • 2
  • 7
  • 15
  • The set of starting values that converge to a particular root is called a basin of attraction. If you google that you'll discover that the maths is very advanced. – Paul Boddington Oct 28 '15 at 22:30
  • 1
    [Are you surprised?](https://en.wikipedia.org/wiki/Newton_fractal) – geometrian Oct 28 '15 at 23:24
  • @PaulBoddington It appears that the OP had actually a bug (or more) in their code rather than a starting value outside the basin of attraction. However, it appears that the OP assumes that Raphson's method (falsely called Newton's -- who failed to publish it) should always converge. – Walter Oct 28 '15 at 23:24

1 Answers1

2

I'm not sure why would you say it's divergent.

I implemented Newton's method similarly to yours:

double g(int w[], int n, double x) {
    double result = 0;
    for (int i = 0; i < n + 1; i++)
        result += w[i] * pow(x, n - i);
    return result;
}

double dg_dx(int w[], int n, double x) {
    double result = 0;
    for (int i = 0; i < n ; i++)
        result += w[i] * (n - i) * pow(x, n - i - 1);
    return result;
}

int main() {

    double xn = 0;        // Choose initial value. I chose 0.
    double gx;
    double dg_dx_x;
    int w[] = { 1, 2, 2, 2, 1 };
    int i = 0;
    int n = 4;

    do {
        gx = g(w, n, xn);
        dg_dx_x = dg_dx(w, n, xn);
        xn = xn - (gx / dg_dx_x);
        i++;
    } while (abs(gx) > 10e-5 && i < 100);

    std::cout << xn << '\n';
}

And it yields -0.997576, which is close to the solution -1.

dorverbin
  • 467
  • 1
  • 4
  • 17
  • hmm, my programm should be flexible for any polynomial, I wrote g and dg, so maybe it causes the problem. Could you look at it? – adolzi Oct 28 '15 at 22:28
  • What did you mean by the variables `w` and `n` in your functions `g` and `dg`? What is the full definition of `g` and `dg`? This should not be a problem if `g` and `dg` are correctly defined, therefore you should post your implementation and this way we could see what went wrong. – dorverbin Oct 28 '15 at 22:32
  • w is the array of coefficients of polynomial and n is the degree. I've just updated my post of dg definition – adolzi Oct 28 '15 at 22:34
  • I wonder if I rather should have used Horner's to compute g and dg ? – adolzi Oct 28 '15 at 22:37
  • I have edited my answer to make the program fit all polynomials using w[] and n. It works fine. The main difference, and I guess this was your issue, is that in the derivative function, which I called `dg_dx`, the stopping condition in the `for` loop should be `i < n`, since the derivative of a constant is zero. – dorverbin Oct 28 '15 at 22:47
  • Actually... Thank you so much – adolzi Oct 28 '15 at 22:54
  • That depends on how you encoded the polynomial. Hopefully with a `double` array for `w` and `n=sizeof(w)/sizeof(w[0]) - 1`. Then the code above converges directly to -1. – Lutz Lehmann Nov 02 '15 at 10:57