Numeric root-finding in Python

[Steven D'Aprano]

This is only peripherally a Python problem, but in case anyone has any 
good ideas I'm going to ask it.

I have a routine to calculate an approximation of Lambert's W function, 
and then apply a root-finding technique to improve the approximation. 
This mostly works well, but sometimes the root-finder gets stuck in a 
cycle.

Here's my function:

import math
def improve(x, w, exp=math.exp):
    """Use Halley's method to improve an estimate of W(x) given 
    an initial estimate w.
    """
    try:
        for i in range(36):  # Max number of iterations.
            ew = exp(w)
            a = w*ew - x
            b = ew*(w + 1)
            err = -a/b  # Estimate of the error in the current w.
            if abs(err) <= 1e-16:
                break
            print '%d: w= %r err= %r' % (i, w, err)
            # Make a better estimate.
            c = (w + 2)*a/(2*w + 2)
            delta = a/(b - c)
            w -= delta
        else:
            raise RuntimeError('calculation failed to converge', err)
    except ZeroDivisionError:
        assert w == -1
    return w


Here's an example where improve() converges very quickly:

py> improve(-0.36, -1.222769842388856)
0: w= -1.222769842388856 err= -2.9158979924038895e-07
1: w= -1.2227701339785069 err= 8.4638038491998997e-16
-1.222770133978506

That's what I expect: convergence in only a few iterations.

Here's an example where it gets stuck in a cycle, bouncing back and forth
between two values:

py> improve(-0.36787344117144249, -1.0057222396915309)
0: w= -1.0057222396915309 err= 2.6521238905750239e-14
1: w= -1.0057222396915044 err= -2.6521238905872001e-14
2: w= -1.0057222396915309 err= 2.6521238905750239e-14
3: w= -1.0057222396915044 err= -2.6521238905872001e-14
4: w= -1.0057222396915309 err= 2.6521238905750239e-14
5: w= -1.0057222396915044 err= -2.6521238905872001e-14
[...]
32: w= -1.0057222396915309 err= 2.6521238905750239e-14
33: w= -1.0057222396915044 err= -2.6521238905872001e-14
34: w= -1.0057222396915309 err= 2.6521238905750239e-14
35: w= -1.0057222396915044 err= -2.6521238905872001e-14
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "<stdin>", line 19, in improve
RuntimeError: ('calculation failed to converge', -2.6521238905872001e-14)

(The correct value for w is approximately -1.00572223991.)

I know that Newton's method is subject to cycles, but I haven't found any 
discussion about Halley's method and cycles, nor do I know what the best 
approach for breaking them would be. None of the papers on calculating
the Lambert W function that I have found mentions this.

Does anyone have any advice for solving this?



-- 
Steven