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.

Are you sure 2000 iterations are enough? From my experience, doing this kind of iterative calculations is better with `while` loops. — Alex Dubrovsky, Feb 19 '18 at 13:43
You should fix the indentation in your code. Python code with obvious indentation errors is not a [mcve]. — khelwood, Feb 19 '18 at 13:44
Fixed those errors, but I've tried a while loop. It does nothing, still stops where it does now, with 4.245665128208564e-05 being the last outputted Er value. — Zach, Feb 19 '18 at 13:51

score 0 · Answer 1 · answered Feb 19 '18 at 13:52

As @alex-dubrovsky metioned, calculation that imply convergence should be implemented with conditional cycles, i.e.:

while True:
    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

With this approach you're always at risk of infinite loop, but this is more or less fine, as algo implies "running until converged"

I just tried that, and I get exactly the same output: 4.245665128208564e-05 for Er and x is 1998.479871524306 — Zach, Feb 19 '18 at 13:55

Alex Dubrovsky · Accepted Answer · 2018-02-19T14:02:20.107

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)

I'll change to while. But I think the professor input the wrong result into the auto-checker, as it checks against the number 1998.47997156. Which is not correct for that function with at least 1.0e-06 accuracy I think — Zach, Feb 19 '18 at 14:01
If you have an access to something like Maple or Mathematica, you can try to solve this equation there and make sure the result is the same as here. — Alex Dubrovsky, Feb 19 '18 at 14:05

Python newton raphson; precision in calculations

2 Answers2