zingy
zingy

Reputation: 811

iterative Newton's method

I have got this code to solve Newton's method for a given polynomial and initial guess value. I want to turn into an iterative process which Newton's method actually is. The program should keeping running till the output value "x_n" becomes constant. And that final value of x_n is the actual root. Also, while using this method in my algorithm it should always produce a positive root between 0 and 1. So does converting the negative output (root) into a positive number would make any difference? Thank you.

import copy

poly = [[-0.25,3], [0.375,2], [-0.375,1], [-3.1,0]]

def poly_diff(poly):
    """ Differentiate a polynomial. """

    newlist = copy.deepcopy(poly)

    for term in newlist:
        term[0] *= term[1]
        term[1] -= 1

    return newlist

def poly_apply(poly, x):
    """ Apply a value to a polynomial. """

    sum = 0.0 

    for term in poly:
        sum += term[0] * (x ** term[1])

    return sum

def poly_root(poly):
    """ Returns a root of the polynomial"""

    poly_d = poly_diff(poly)
    x = float(raw_input("Enter initial guess:"))

    x_n = x - (float(poly_apply(poly, x)) / poly_apply(poly_d, x))

    print x_n

if __name__ == "__main__" :
    poly_root(poly)

Upvotes: 1

Views: 1788

Answers (2)

marq
marq

Reputation: 1

import copy

poly = [[1,64], [2,109], [3,137], [4,138], [5,171], [6,170]]

def poly_diff(poly):

    newlist = copy.deepcopy(poly)

    for term in newlist:
        term[0] *= term[1]
        term[1] -= 1

    return newlist

def poly_apply(poly, x):

    sum = 0.0 

    for term in poly:
        sum += term[0] * (x ** term[1])

    return sum

def poly_root(poly):

    poly_d = poly_diff(poly)
    x = float(input("Enter initial guess:"))

    x_n = x - (float(poly_apply(poly, x)) / poly_apply(poly_d, x))

    print (x_n)

if __name__ == "__main__" :
    poly_root(poly)

Upvotes: 0

wberry
wberry

Reputation: 19337

First, in poly_diff, you should check to see if the exponent is zero, and if so simply remove that term from the result. Otherwise you will end up with the derivative being undefined at zero.

def poly_root(poly):
    """ Returns a root of the polynomial"""
    poly_d = poly_diff(poly)
    x = None
    x_n = float(raw_input("Enter initial guess:"))
    while x != x_n:
        x = x_n
        x_n = x - (float(poly_apply(poly, x)) / poly_apply(poly_d, x))
    return x_n

That should do it. However, I think it is possible that for certain polynomials this may not terminate, due to floating point rounding error. It may end up in a repeating cycle of approximations that differ only in the least significant bits. You might terminate when the percentage of change reaches a lower limit, or after a number of iterations.

Upvotes: 1

Related Questions