Welcome, guest | Sign In | My Account | Store | Cart

Halley's method for solving equations of type f(x)=0.

Python, 29 lines
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
# 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

def fp(x):
    global h
    return (f(x + h) - f(x)) / h

def fpp(x):
    global h
    return (fp(x + h) - fp(x)) / h

# main
x = 2.0 # initial value

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

2 comments

FB36 (author) 11 years ago  # | flag

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
FB36 (author) 11 years ago  # | flag

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)
Created by FB36 on Sat, 27 Nov 2010 (MIT)
Python recipes (4591)
FB36's recipes (148)

Required Modules

  • (none specified)

Other Information and Tasks