residues.py 2.2 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273
  1. """
  2. This module implements the Residue function and related tools for working
  3. with residues.
  4. """
  5. from sympy.core.mul import Mul
  6. from sympy.core.singleton import S
  7. from sympy.core.sympify import sympify
  8. from sympy.utilities.timeutils import timethis
  9. @timethis('residue')
  10. def residue(expr, x, x0):
  11. """
  12. Finds the residue of ``expr`` at the point x=x0.
  13. The residue is defined as the coefficient of ``1/(x-x0)`` in the power series
  14. expansion about ``x=x0``.
  15. Examples
  16. ========
  17. >>> from sympy import Symbol, residue, sin
  18. >>> x = Symbol("x")
  19. >>> residue(1/x, x, 0)
  20. 1
  21. >>> residue(1/x**2, x, 0)
  22. 0
  23. >>> residue(2/sin(x), x, 0)
  24. 2
  25. This function is essential for the Residue Theorem [1].
  26. References
  27. ==========
  28. .. [1] https://en.wikipedia.org/wiki/Residue_theorem
  29. """
  30. # The current implementation uses series expansion to
  31. # calculate it. A more general implementation is explained in
  32. # the section 5.6 of the Bronstein's book {M. Bronstein:
  33. # Symbolic Integration I, Springer Verlag (2005)}. For purely
  34. # rational functions, the algorithm is much easier. See
  35. # sections 2.4, 2.5, and 2.7 (this section actually gives an
  36. # algorithm for computing any Laurent series coefficient for
  37. # a rational function). The theory in section 2.4 will help to
  38. # understand why the resultant works in the general algorithm.
  39. # For the definition of a resultant, see section 1.4 (and any
  40. # previous sections for more review).
  41. from sympy.series.order import Order
  42. from sympy.simplify.radsimp import collect
  43. expr = sympify(expr)
  44. if x0 != 0:
  45. expr = expr.subs(x, x + x0)
  46. for n in (0, 1, 2, 4, 8, 16, 32):
  47. s = expr.nseries(x, n=n)
  48. if not s.has(Order) or s.getn() >= 0:
  49. break
  50. s = collect(s.removeO(), x)
  51. if s.is_Add:
  52. args = s.args
  53. else:
  54. args = [s]
  55. res = S.Zero
  56. for arg in args:
  57. c, m = arg.as_coeff_mul(x)
  58. m = Mul(*m)
  59. if not (m in (S.One, x) or (m.is_Pow and m.exp.is_Integer)):
  60. raise NotImplementedError('term of unexpected form: %s' % m)
  61. if m == 1/x:
  62. res += c
  63. return res