trigonometry.py 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335
  1. from sympy.core import cacheit, Dummy, Ne, Integer, Rational, S, Wild
  2. from sympy.functions import binomial, sin, cos, Piecewise
  3. from .integrals import integrate
  4. # TODO sin(a*x)*cos(b*x) -> sin((a+b)x) + sin((a-b)x) ?
  5. # creating, each time, Wild's and sin/cos/Mul is expensive. Also, our match &
  6. # subs are very slow when not cached, and if we create Wild each time, we
  7. # effectively block caching.
  8. #
  9. # so we cache the pattern
  10. # need to use a function instead of lamda since hash of lambda changes on
  11. # each call to _pat_sincos
  12. def _integer_instance(n):
  13. return isinstance(n, Integer)
  14. @cacheit
  15. def _pat_sincos(x):
  16. a = Wild('a', exclude=[x])
  17. n, m = [Wild(s, exclude=[x], properties=[_integer_instance])
  18. for s in 'nm']
  19. pat = sin(a*x)**n * cos(a*x)**m
  20. return pat, a, n, m
  21. _u = Dummy('u')
  22. def trigintegrate(f, x, conds='piecewise'):
  23. """
  24. Integrate f = Mul(trig) over x.
  25. Examples
  26. ========
  27. >>> from sympy import sin, cos, tan, sec
  28. >>> from sympy.integrals.trigonometry import trigintegrate
  29. >>> from sympy.abc import x
  30. >>> trigintegrate(sin(x)*cos(x), x)
  31. sin(x)**2/2
  32. >>> trigintegrate(sin(x)**2, x)
  33. x/2 - sin(x)*cos(x)/2
  34. >>> trigintegrate(tan(x)*sec(x), x)
  35. 1/cos(x)
  36. >>> trigintegrate(sin(x)*tan(x), x)
  37. -log(sin(x) - 1)/2 + log(sin(x) + 1)/2 - sin(x)
  38. References
  39. ==========
  40. .. [1] http://en.wikibooks.org/wiki/Calculus/Integration_techniques
  41. See Also
  42. ========
  43. sympy.integrals.integrals.Integral.doit
  44. sympy.integrals.integrals.Integral
  45. """
  46. pat, a, n, m = _pat_sincos(x)
  47. f = f.rewrite('sincos')
  48. M = f.match(pat)
  49. if M is None:
  50. return
  51. n, m = M[n], M[m]
  52. if n.is_zero and m.is_zero:
  53. return x
  54. zz = x if n.is_zero else S.Zero
  55. a = M[a]
  56. if n.is_odd or m.is_odd:
  57. u = _u
  58. n_, m_ = n.is_odd, m.is_odd
  59. # take smallest n or m -- to choose simplest substitution
  60. if n_ and m_:
  61. # Make sure to choose the positive one
  62. # otherwise an incorrect integral can occur.
  63. if n < 0 and m > 0:
  64. m_ = True
  65. n_ = False
  66. elif m < 0 and n > 0:
  67. n_ = True
  68. m_ = False
  69. # Both are negative so choose the smallest n or m
  70. # in absolute value for simplest substitution.
  71. elif (n < 0 and m < 0):
  72. n_ = n > m
  73. m_ = not (n > m)
  74. # Both n and m are odd and positive
  75. else:
  76. n_ = (n < m) # NB: careful here, one of the
  77. m_ = not (n < m) # conditions *must* be true
  78. # n m u=C (n-1)/2 m
  79. # S(x) * C(x) dx --> -(1-u^2) * u du
  80. if n_:
  81. ff = -(1 - u**2)**((n - 1)/2) * u**m
  82. uu = cos(a*x)
  83. # n m u=S n (m-1)/2
  84. # S(x) * C(x) dx --> u * (1-u^2) du
  85. elif m_:
  86. ff = u**n * (1 - u**2)**((m - 1)/2)
  87. uu = sin(a*x)
  88. fi = integrate(ff, u) # XXX cyclic deps
  89. fx = fi.subs(u, uu)
  90. if conds == 'piecewise':
  91. return Piecewise((fx / a, Ne(a, 0)), (zz, True))
  92. return fx / a
  93. # n & m are both even
  94. #
  95. # 2k 2m 2l 2l
  96. # we transform S (x) * C (x) into terms with only S (x) or C (x)
  97. #
  98. # example:
  99. # 100 4 100 2 2 100 4 2
  100. # S (x) * C (x) = S (x) * (1-S (x)) = S (x) * (1 + S (x) - 2*S (x))
  101. #
  102. # 104 102 100
  103. # = S (x) - 2*S (x) + S (x)
  104. # 2k
  105. # then S is integrated with recursive formula
  106. # take largest n or m -- to choose simplest substitution
  107. n_ = (abs(n) > abs(m))
  108. m_ = (abs(m) > abs(n))
  109. res = S.Zero
  110. if n_:
  111. # 2k 2 k i 2i
  112. # C = (1 - S ) = sum(i, (-) * B(k, i) * S )
  113. if m > 0:
  114. for i in range(0, m//2 + 1):
  115. res += (S.NegativeOne**i * binomial(m//2, i) *
  116. _sin_pow_integrate(n + 2*i, x))
  117. elif m == 0:
  118. res = _sin_pow_integrate(n, x)
  119. else:
  120. # m < 0 , |n| > |m|
  121. # /
  122. # |
  123. # | m n
  124. # | cos (x) sin (x) dx =
  125. # |
  126. # |
  127. #/
  128. # /
  129. # |
  130. # -1 m+1 n-1 n - 1 | m+2 n-2
  131. # ________ cos (x) sin (x) + _______ | cos (x) sin (x) dx
  132. # |
  133. # m + 1 m + 1 |
  134. # /
  135. res = (Rational(-1, m + 1) * cos(x)**(m + 1) * sin(x)**(n - 1) +
  136. Rational(n - 1, m + 1) *
  137. trigintegrate(cos(x)**(m + 2)*sin(x)**(n - 2), x))
  138. elif m_:
  139. # 2k 2 k i 2i
  140. # S = (1 - C ) = sum(i, (-) * B(k, i) * C )
  141. if n > 0:
  142. # / /
  143. # | |
  144. # | m n | -m n
  145. # | cos (x)*sin (x) dx or | cos (x) * sin (x) dx
  146. # | |
  147. # / /
  148. #
  149. # |m| > |n| ; m, n >0 ; m, n belong to Z - {0}
  150. # n 2
  151. # sin (x) term is expanded here in terms of cos (x),
  152. # and then integrated.
  153. #
  154. for i in range(0, n//2 + 1):
  155. res += (S.NegativeOne**i * binomial(n//2, i) *
  156. _cos_pow_integrate(m + 2*i, x))
  157. elif n == 0:
  158. # /
  159. # |
  160. # | 1
  161. # | _ _ _
  162. # | m
  163. # | cos (x)
  164. # /
  165. #
  166. res = _cos_pow_integrate(m, x)
  167. else:
  168. # n < 0 , |m| > |n|
  169. # /
  170. # |
  171. # | m n
  172. # | cos (x) sin (x) dx =
  173. # |
  174. # |
  175. #/
  176. # /
  177. # |
  178. # 1 m-1 n+1 m - 1 | m-2 n+2
  179. # _______ cos (x) sin (x) + _______ | cos (x) sin (x) dx
  180. # |
  181. # n + 1 n + 1 |
  182. # /
  183. res = (Rational(1, n + 1) * cos(x)**(m - 1)*sin(x)**(n + 1) +
  184. Rational(m - 1, n + 1) *
  185. trigintegrate(cos(x)**(m - 2)*sin(x)**(n + 2), x))
  186. else:
  187. if m == n:
  188. ##Substitute sin(2x)/2 for sin(x)cos(x) and then Integrate.
  189. res = integrate((sin(2*x)*S.Half)**m, x)
  190. elif (m == -n):
  191. if n < 0:
  192. # Same as the scheme described above.
  193. # the function argument to integrate in the end will
  194. # be 1, this cannot be integrated by trigintegrate.
  195. # Hence use sympy.integrals.integrate.
  196. res = (Rational(1, n + 1) * cos(x)**(m - 1) * sin(x)**(n + 1) +
  197. Rational(m - 1, n + 1) *
  198. integrate(cos(x)**(m - 2) * sin(x)**(n + 2), x))
  199. else:
  200. res = (Rational(-1, m + 1) * cos(x)**(m + 1) * sin(x)**(n - 1) +
  201. Rational(n - 1, m + 1) *
  202. integrate(cos(x)**(m + 2)*sin(x)**(n - 2), x))
  203. if conds == 'piecewise':
  204. return Piecewise((res.subs(x, a*x) / a, Ne(a, 0)), (zz, True))
  205. return res.subs(x, a*x) / a
  206. def _sin_pow_integrate(n, x):
  207. if n > 0:
  208. if n == 1:
  209. #Recursion break
  210. return -cos(x)
  211. # n > 0
  212. # / /
  213. # | |
  214. # | n -1 n-1 n - 1 | n-2
  215. # | sin (x) dx = ______ cos (x) sin (x) + _______ | sin (x) dx
  216. # | |
  217. # | n n |
  218. #/ /
  219. #
  220. #
  221. return (Rational(-1, n) * cos(x) * sin(x)**(n - 1) +
  222. Rational(n - 1, n) * _sin_pow_integrate(n - 2, x))
  223. if n < 0:
  224. if n == -1:
  225. ##Make sure this does not come back here again.
  226. ##Recursion breaks here or at n==0.
  227. return trigintegrate(1/sin(x), x)
  228. # n < 0
  229. # / /
  230. # | |
  231. # | n 1 n+1 n + 2 | n+2
  232. # | sin (x) dx = _______ cos (x) sin (x) + _______ | sin (x) dx
  233. # | |
  234. # | n + 1 n + 1 |
  235. #/ /
  236. #
  237. return (Rational(1, n + 1) * cos(x) * sin(x)**(n + 1) +
  238. Rational(n + 2, n + 1) * _sin_pow_integrate(n + 2, x))
  239. else:
  240. #n == 0
  241. #Recursion break.
  242. return x
  243. def _cos_pow_integrate(n, x):
  244. if n > 0:
  245. if n == 1:
  246. #Recursion break.
  247. return sin(x)
  248. # n > 0
  249. # / /
  250. # | |
  251. # | n 1 n-1 n - 1 | n-2
  252. # | sin (x) dx = ______ sin (x) cos (x) + _______ | cos (x) dx
  253. # | |
  254. # | n n |
  255. #/ /
  256. #
  257. return (Rational(1, n) * sin(x) * cos(x)**(n - 1) +
  258. Rational(n - 1, n) * _cos_pow_integrate(n - 2, x))
  259. if n < 0:
  260. if n == -1:
  261. ##Recursion break
  262. return trigintegrate(1/cos(x), x)
  263. # n < 0
  264. # / /
  265. # | |
  266. # | n -1 n+1 n + 2 | n+2
  267. # | cos (x) dx = _______ sin (x) cos (x) + _______ | cos (x) dx
  268. # | |
  269. # | n + 1 n + 1 |
  270. #/ /
  271. #
  272. return (Rational(-1, n + 1) * sin(x) * cos(x)**(n + 1) +
  273. Rational(n + 2, n + 1) * _cos_pow_integrate(n + 2, x))
  274. else:
  275. # n == 0
  276. #Recursion Break.
  277. return x