Halley's Method for Solving Equations (Python recipe) by FB36

ActiveState Code (http://code.activestate.com/recipes/577472/)

Here is another version that uses a generalized derivative function:

# Halley's method for solving f(x)=0
# http://en.wikipedia.org/wiki/Halley%27s_method
# FB - 201011265
global h
h = 0.00000001
eps = 0.000001

# f(x) to solve
def f(x):
    return x * x - 2.0

# general numerical derivative function
def fp(x, k):
    global h
    if k == 0:
        return f(x)
    else:
        return (fp(x + h, k - 1) - fp(x, k - 1)) / h

# main
x = 2.0 # initial value

while True:
    fx = f(x)
    fpx = fp(x, 1)
    xnew = x - (2.0 * fx * fpx) / (2.0 * fpx * fpx - fx * fp(x, 2))
    print xnew
    if abs(xnew - x) <= eps: break
    x = xnew

Numerical derivative function can be improved even further:

# generalized numerical derivative function
def fp(x, k):
    global h
    if k == 0:
        return f(x)
    else:
        # two-point estimation
        # return (fp(x + h, k - 1) - fp(x, k - 1)) / h
        # three-point estimation
        # return (fp(x + h, k - 1) - fp(x - h, k - 1)) / (2 * h)
        # five-point stencil
        return (-fp(x+2*h,k-1)+8*fp(x+h,k-1)-8*fp(x-h,k-1)+fp(x-2*h,k-1))/(12*h)