◐
Written by Niels Moseley
on 

Newton's and Halley's methods

The following Python3 code shows how to use Newton and Halley’s methods to iteratively calculate the $x$ for which a function $f(x)$ is zero. This can be used to calculate square roots and other things, such as $\pi$.

The function $f(x)$ must be continuous (the first and second order derivatives must exist) over the search area for these methods to work. Also, the initial starting point / guess for $x$ should be sufficiently close to the desired solution.


import math

# Numerically calculate the derivative of a function fx
# at point x
def deriv(fx, x : float) -> float:
    dx = 0.00001
    return 0.5 * (fx(x+dx)-fx(x-dx)) / dx

# Numerically calculate the second derivative of a function fx
# at point x
def deriv2(fx, x : float) -> float:
    dx = 0.00001
    return 0.5 * (deriv(fx, x+dx) -deriv(fx, x-dx)) / dx
    
# Newton's method:
# x(n+1) = x(n) - fn(x(n)) / fn'(x(n))
# f(x) = x^2 - c
# 
# analytical version for f(x) = x^2 - c
# -> f'(x) = 2x
# so newton update rule is
# x(n+1) = x(n) - (x(n)^2 - c) / 2x(n)
# Note: the above derivation is not used here.
#
def newton(fx, guess : float) -> float:
    print("  guess:", guess)
    newguess = guess - fx(guess) / deriv(fx, guess)
    tol = 1.0e-12
    err = fx(newguess)
    if (abs(err) < tol):
        return newguess;
    else:
        return newton(fx, newguess)
        
# Halley's method:
# x(n+1) = x(n) - 2[fn(x(n))*fn'(x(n))] / ( 2fn'(x(n))^2 - fn(x(n))fn''(x(n))
# f(x) = x^2 - c
# 
# analytical version for f(x) = x^2 - c
# -> f'(x) = 2x
# -> f''(x) = 2
# so halley update rule is
# x(n+1) = x(n) - 2[(x^2 - c)*(2x)]/( 2(2x)^2 - (x^2-c)*2)
# x(n+1) = x(n) - (4x^3 - 4cx)/(8x^2 - 2x^2 -2c)
# Note: the above derivation is not used here.
#
def halley(fx, guess : float) -> float:
    v =  fx(guess)
    d  = deriv(fx, guess)
    d2 = deriv2(fx, guess)
    print("  guess:", guess, "  d2:", d2)
    newguess = guess - (2.0*v*d) / (2.0*d*d - v*d2)
    tol = 1.0e-10
    err = fx(newguess)
    if (abs(err) < tol):
        return newguess;
    else:
        return halley(fx, newguess)

# Function is zero when x equals sqrt(2)
def func(x : float) -> float:
    return x*x - 2.0

# Function is zero when x equals pi
def func2(x : float) -> float:
    return math.sin(x)

# Function is zero when x equals pi/4
def func3(x : float) -> float:
    return math.tan(x)-1.0

# Use Netwon and Halley to calculate sqrt(2)
print(newton(func, 1.0))
print(halley(func, 1.0))
print("Ref: 1.41421356237")

# Use Netwon and Halley to calculate pi
print(newton(func2, 2))
print(halley(func2, 2))
print("Ref: 3.14159265359")

# Use Netwon and Halley to calculate pi/4
print(newton(func3, 0.75)*4)
print(halley(func3, 0.75)*4)
print("Ref: 3.14159265359")

Output


  guess: 1.0
  guess: 1.4999999999995
  guess: 1.4166666666672125
  guess: 1.4142156862745259
  guess: 1.4142135623746899
1.4142135623730951
  guess: 1.0   d2: 2.000001275714869
  guess: 1.3999999489711716   d2: 2.000000165480742
  guess: 1.4142131979595292   d2: 2.0000007205922543
1.4142135623730951
Ref:  1.41421356237
  guess: 2
  guess: 4.1850398632966375
  guess: 2.4678936743166275
  guess: 3.2661862777023662
  guess: 3.1409439123142335
  guess: 3.1415926536808105
3.141592653589793
  guess: 2   d2: -0.909297359630301
  guess: 2.645087456831022   d2: -0.47635576039262156
  guess: 3.11752605316398   d2: -0.02406428355472556
  guess: 3.1415903294119794   d2: -2.3241797375561646e-06
3.141592653589793
Ref:  3.141592653589793238
  guess: 0.75
  guess: 0.7866211075276159
  guess: 0.7853996602085616
  guess: 0.7853981633996889
3.141592653589793
  guess: 0.75   d2: 3.4802058457295
  guess: 0.7854129557853823   d2: 4.000236530909973
3.1415926535897962
Ref:  3.141592653589793238

Further reading