Newton Method: building higher order procedure using python

Question

I am working through Structure and Interpretation of Computer Programs.

in pg. 73, it uses Newton's Method as example of how to construct higher order procedure.

here is my code:

def deriv(g):
    dx = 0.00001
    return lambda x: (g(x + dx) - g(x)) / dx


def newton_transform(g):
    return lambda x: x - g(x) / deriv(g)(x)


def fixed_point(f, guess):
    def close_enough(a, b):
        tolerance = 0.00001
        return abs(a - b) < tolerance

    def a_try(guess):
        next = f(guess)
        if close_enough(guess, next):
            return next
        else:
            return a_try(next)

    return a_try(guess)


def newton_method(g, guess):
    return fixed_point(newton_transform(g), guess)


def sqrt(x):
    return newton_method(lambda y: x / y, 1.0)

print sqrt(2)

The code will crash and give me ZeroDivisionError. I know how it crash but I don't understand why it behave as it is.

In my "a_try" function, every "next" is doubled value of "guess". When I print out "guess" and "next" of every iteration, my next guess just keep doubling. So the integer overflow in the end.

Why? What's wrong with my code? What's wrong with my logic? Thanks for your time. Please help.

Answer 1

To use Newton's method to find, say, sqrt(2)—that is, y**2 == 2—you first write g(y) = y**2 - 2, and then iterate over that with newton_transform until it converges.

Your deriv and newton_transform are fine, and your fixed_point does in fact iterate over newton_transform until it converges—or until you hit the recursion limit, or until you underflow a float operation. In your case, it's the last one.

Why? Well, look at your g(y): it's 2/y. I don't know where you got that from, but g / g' is just -y, so the newton transform is y - -y, which obviously isn't going to converge.

But if you plug in y**2 - 2, then the transform on g/g' will converge (at least for most values).

So:

def sqrt(x):
    return newton_method(lambda y: y**2-x, 1.0)