Python newton raphson; precision in calculations

Question

I've written this function in python:

def f2(x):
    return (5.0*x + log1p(x) - 10000.0)

def dfdx2(x):
    return (5.0-(1.0/x))

def newtonRaphson2(f, dfdx, x, tol):

    x0 = x

    for i in range(1, 2000):
        if f(x) == 0.0:
            return x

    if dfdx(x) == 0.0:
        print(dfdx(x))
        break

    x = x - (f(x) / dfdx(x))
    #print(x)

    Er = abs(x0-x)/abs(x0)
    if Er <= tol:
        return x
    print(Er)
    x0 = x
    return x

Then I execute it like this:

task2 = newtonRaphson2(f2, dfdx2, 1, 0.000001);
print(task2)

For the output check the Er outputs final accuracy of 4.245665128208564e-05 before it returns x.

X is returned at 1998.479871524306, which is a pretty good estimate, but I'm preferably looking to get it down to 1.0e-06 at least. Changing tol variable to 1.0e-08 seems to do nothing.

I'm guessing maybe putting every variable into double is a better idea, but I still have no idea why my code stops where it does. I'm not that stable with python either, which is why I'm asking. I've already written one of these that works, but its for a far simpler equation.

Answer 1

Your code works fine, once you indent it properly and add from math import log1p . Just put the print(Er) line right after it's calculated to see its final value. Er gets to ~10^-9. This worked for me:

from math import log1p

def f2(x):
    return (5.0*x + log1p(x) - 10000.0)

def dfdx2(x):
    return (5.0-(1.0/x))

def newtonRaphson2(f, dfdx, x, tol):

    x0 = x
    for i in range(1, 2000):
        if f(x) == 0.0:
            return x

        if dfdx(x) == 0.0:
            print(dfdx(x))
            break
        x = x - (f(x) / dfdx(x))
        #print(x)

        Er = abs(x0-x)/abs(x0)
        print('Er = {}'.format(Er))
        if Er <= tol:
            return x

        x0 = x
    return x

x = newtonRaphson2(f2, dfdx2, 1, 0.0000001)
print 'X = {}'.format(x)

The output was:

Er = 2498.5767132
Er = 0.200506616666
Er = 4.24566512821e-05
Er = 8.49642413214e-09
X = 1998.47987152

Consider using while here. Newton-Raphson algorithm typically converges very fast, so you won't need many iterations to run most cases.

This gives the same result:

    from math import log1p

def f2(x):
    return (5.0*x + log1p(x) - 10000.0)

def dfdx2(x):
    return (5.0-(1.0/x))

def newtonRaphson2(f, dfdx, x, tol):

    x0 = x
    Er = 1
    while Er >= tol:
        if f(x) == 0.0:
            return x

        if dfdx(x) == 0.0:
            print(dfdx(x))
            break
        x = x - (f(x) / dfdx(x))
        #print(x)

        Er = abs(x0-x)/abs(x0)
        print('Er = {}'.format(Er))

        x0 = x
    return x

x = newtonRaphson2(f2, dfdx2, 1, 0.0000001)
print 'X = {}'.format(x)