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.

Integers don't overflow in Python. They just keep growing. (Well, eventually you get a `MemoryError`, but they don't roll over to 0.) So, your description of the problem can't be right. However, you're obviously converting to float somewhere (hopefully early enough) or `sqrt(2)` would have to be `1`, so this isn't really relevant. — abarnert, Nov 14 '14 at 01:05
Maybe add `from __future__ import division` at the top work? — wilbeibi, Nov 14 '14 at 01:08
I'm confused about how this is supposed to work. There's no squaring happening anywhere, so, even if this did converge, how would it converge to the square root? (Your `deriv` and `newton_transform` functions look fine, I just don't get what you're trying to transform here.) Without seeing the Scheme code you're porting it's hard to guess what's missing — abarnert, Nov 14 '14 at 01:18

abarnert · Accepted Answer · 2014-11-14T01:37:44.207

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)

edited Nov 14 '14 at 01:37

answered Nov 14 '14 at 01:30

abarnert

354,177
51
601
671

Yes. This is the answer. I understand newton's method wrong. Thank you very much !!!!!! Logical bug is so much harder to identify than programming bug... – Lucas Shen Nov 14 '14 at 01:33
@LucasShen: I'm glad you understood, because I'm having a really hard time trying to explain the math, and I was going to suggest that if you didn't get it you should go somewhere else (maybe [Math.sx](http://mathematics.stackexchange.com)) where someone can put it better. :) – abarnert Nov 14 '14 at 01:35
hahaha, that's good resource too. But just to your reference, from my previous mistake, my g(y) = x/y, so my newton's transform will get me a lambda y: 2*y. That's the reason my guess couldn't converge...and keep doubling every iteration. – Lucas Shen Nov 14 '14 at 01:37

Newton Method: building higher order procedure using python

1 Answers1