Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.python > #20275

Re: Numeric root-finding in Python

Show key headers only | View raw

I don't know the first thing about this math problem however,

if I were to code this I might try ;

   except ZeroDivisionError:
         assert w = -1

rather than;

   except ZeroDivisionError:
         assert w == -1


jimonlinux



On Sunday, February 12, 2012 06:41:20 AM Steven D'Aprano wrote:
> 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?

Back to comp.lang.python | Previous | Next — Previous in thread | Next in thread | Find similar | Unroll thread

Thread

Numeric root-finding in Python Steven D'Aprano <steve+comp.lang.python@pearwood.info> - 2012-02-12 06:41 +0000
  Re: Numeric root-finding in Python Eelco <hoogendoorn.eelco@gmail.com> - 2012-02-12 02:10 -0800
    Re: Numeric root-finding in Python Mark Lawrence <breamoreboy@yahoo.co.uk> - 2012-02-12 13:39 +0000
    Re: Numeric root-finding in Python Terry Reedy <tjreedy@udel.edu> - 2012-02-12 16:19 -0500
  Re: Numeric root-finding in Python 88888 Dihedral <dihedral88888@googlemail.com> - 2012-02-12 05:13 -0800
  Re: Numeric root-finding in Python Robert Kern <robert.kern@gmail.com> - 2012-02-12 13:52 +0000
    Re: Numeric root-finding in Python Steven D'Aprano <steve+comp.lang.python@pearwood.info> - 2012-02-12 22:53 +0000
  Re: Numeric root-finding in Python inq1ltd <inq1ltd@inqvista.com> - 2012-02-12 10:20 -0500
    Re: Numeric root-finding in Python David Monaghan <monaghand.david@gmail.com> - 2012-02-21 01:33 +0000
      Re: Numeric root-finding in Python Mark Lawrence <breamoreboy@yahoo.co.uk> - 2012-02-21 02:21 +0000
  Re: Numeric root-finding in Python Dennis Lee Bieber <wlfraed@ix.netcom.com> - 2012-02-12 11:47 -0500
  Re: Numeric root-finding in Python Dave Angel <d@davea.name> - 2012-02-12 12:00 -0500
  Re: Numeric root-finding in Python Mark Dickinson <mdickinson@enthought.com> - 2012-02-12 12:18 -0800
    Re: Numeric root-finding in Python Steven D'Aprano <steve+comp.lang.python@pearwood.info> - 2012-02-12 23:05 +0000
      Re: Numeric root-finding in Python Dave Angel <d@davea.name> - 2012-02-12 18:52 -0500

csiph-web