_lambertw.py 3.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149
  1. from ._ufuncs import _lambertw
  2. import numpy as np
  3. def lambertw(z, k=0, tol=1e-8):
  4. r"""
  5. lambertw(z, k=0, tol=1e-8)
  6. Lambert W function.
  7. The Lambert W function `W(z)` is defined as the inverse function
  8. of ``w * exp(w)``. In other words, the value of ``W(z)`` is
  9. such that ``z = W(z) * exp(W(z))`` for any complex number
  10. ``z``.
  11. The Lambert W function is a multivalued function with infinitely
  12. many branches. Each branch gives a separate solution of the
  13. equation ``z = w exp(w)``. Here, the branches are indexed by the
  14. integer `k`.
  15. Parameters
  16. ----------
  17. z : array_like
  18. Input argument.
  19. k : int, optional
  20. Branch index.
  21. tol : float, optional
  22. Evaluation tolerance.
  23. Returns
  24. -------
  25. w : array
  26. `w` will have the same shape as `z`.
  27. See Also
  28. --------
  29. wrightomega : the Wright Omega function
  30. Notes
  31. -----
  32. All branches are supported by `lambertw`:
  33. * ``lambertw(z)`` gives the principal solution (branch 0)
  34. * ``lambertw(z, k)`` gives the solution on branch `k`
  35. The Lambert W function has two partially real branches: the
  36. principal branch (`k = 0`) is real for real ``z > -1/e``, and the
  37. ``k = -1`` branch is real for ``-1/e < z < 0``. All branches except
  38. ``k = 0`` have a logarithmic singularity at ``z = 0``.
  39. **Possible issues**
  40. The evaluation can become inaccurate very close to the branch point
  41. at ``-1/e``. In some corner cases, `lambertw` might currently
  42. fail to converge, or can end up on the wrong branch.
  43. **Algorithm**
  44. Halley's iteration is used to invert ``w * exp(w)``, using a first-order
  45. asymptotic approximation (O(log(w)) or `O(w)`) as the initial estimate.
  46. The definition, implementation and choice of branches is based on [2]_.
  47. References
  48. ----------
  49. .. [1] https://en.wikipedia.org/wiki/Lambert_W_function
  50. .. [2] Corless et al, "On the Lambert W function", Adv. Comp. Math. 5
  51. (1996) 329-359.
  52. https://cs.uwaterloo.ca/research/tr/1993/03/W.pdf
  53. Examples
  54. --------
  55. The Lambert W function is the inverse of ``w exp(w)``:
  56. >>> import numpy as np
  57. >>> from scipy.special import lambertw
  58. >>> w = lambertw(1)
  59. >>> w
  60. (0.56714329040978384+0j)
  61. >>> w * np.exp(w)
  62. (1.0+0j)
  63. Any branch gives a valid inverse:
  64. >>> w = lambertw(1, k=3)
  65. >>> w
  66. (-2.8535817554090377+17.113535539412148j)
  67. >>> w*np.exp(w)
  68. (1.0000000000000002+1.609823385706477e-15j)
  69. **Applications to equation-solving**
  70. The Lambert W function may be used to solve various kinds of
  71. equations. We give two examples here.
  72. First, the function can be used to solve implicit equations of the
  73. form
  74. :math:`x = a + b e^{c x}`
  75. for :math:`x`. We assume :math:`c` is not zero. After a little
  76. algebra, the equation may be written
  77. :math:`z e^z = -b c e^{a c}`
  78. where :math:`z = c (a - x)`. :math:`z` may then be expressed using
  79. the Lambert W function
  80. :math:`z = W(-b c e^{a c})`
  81. giving
  82. :math:`x = a - W(-b c e^{a c})/c`
  83. For example,
  84. >>> a = 3
  85. >>> b = 2
  86. >>> c = -0.5
  87. The solution to :math:`x = a + b e^{c x}` is:
  88. >>> x = a - lambertw(-b*c*np.exp(a*c))/c
  89. >>> x
  90. (3.3707498368978794+0j)
  91. Verify that it solves the equation:
  92. >>> a + b*np.exp(c*x)
  93. (3.37074983689788+0j)
  94. The Lambert W function may also be used find the value of the infinite
  95. power tower :math:`z^{z^{z^{\ldots}}}`:
  96. >>> def tower(z, n):
  97. ... if n == 0:
  98. ... return z
  99. ... return z ** tower(z, n-1)
  100. ...
  101. >>> tower(0.5, 100)
  102. 0.641185744504986
  103. >>> -lambertw(-np.log(0.5)) / np.log(0.5)
  104. (0.64118574450498589+0j)
  105. """
  106. # TODO: special expert should inspect this
  107. # interception; better place to do it?
  108. k = np.asarray(k, dtype=np.dtype("long"))
  109. return _lambertw(z, k, tol)