transforms.py 50 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588
  1. """ Integral Transforms """
  2. from functools import reduce, wraps
  3. from itertools import repeat
  4. from sympy.core import S, pi
  5. from sympy.core.add import Add
  6. from sympy.core.function import (
  7. AppliedUndef, count_ops, expand, expand_mul, Function)
  8. from sympy.core.mul import Mul
  9. from sympy.core.numbers import igcd, ilcm
  10. from sympy.core.sorting import default_sort_key
  11. from sympy.core.symbol import Dummy
  12. from sympy.core.traversal import postorder_traversal
  13. from sympy.functions.combinatorial.factorials import factorial, rf
  14. from sympy.functions.elementary.complexes import re, arg, Abs
  15. from sympy.functions.elementary.exponential import exp, exp_polar
  16. from sympy.functions.elementary.hyperbolic import cosh, coth, sinh, tanh
  17. from sympy.functions.elementary.integers import ceiling
  18. from sympy.functions.elementary.miscellaneous import Max, Min, sqrt
  19. from sympy.functions.elementary.piecewise import piecewise_fold
  20. from sympy.functions.elementary.trigonometric import cos, cot, sin, tan
  21. from sympy.functions.special.bessel import besselj
  22. from sympy.functions.special.delta_functions import Heaviside
  23. from sympy.functions.special.gamma_functions import gamma
  24. from sympy.functions.special.hyper import meijerg
  25. from sympy.integrals import integrate, Integral
  26. from sympy.integrals.meijerint import _dummy
  27. from sympy.logic.boolalg import to_cnf, conjuncts, disjuncts, Or, And
  28. from sympy.polys.polyroots import roots
  29. from sympy.polys.polytools import factor, Poly
  30. from sympy.polys.rootoftools import CRootOf
  31. from sympy.utilities.iterables import iterable
  32. from sympy.utilities.misc import debug
  33. ##########################################################################
  34. # Helpers / Utilities
  35. ##########################################################################
  36. class IntegralTransformError(NotImplementedError):
  37. """
  38. Exception raised in relation to problems computing transforms.
  39. Explanation
  40. ===========
  41. This class is mostly used internally; if integrals cannot be computed
  42. objects representing unevaluated transforms are usually returned.
  43. The hint ``needeval=True`` can be used to disable returning transform
  44. objects, and instead raise this exception if an integral cannot be
  45. computed.
  46. """
  47. def __init__(self, transform, function, msg):
  48. super().__init__(
  49. "%s Transform could not be computed: %s." % (transform, msg))
  50. self.function = function
  51. class IntegralTransform(Function):
  52. """
  53. Base class for integral transforms.
  54. Explanation
  55. ===========
  56. This class represents unevaluated transforms.
  57. To implement a concrete transform, derive from this class and implement
  58. the ``_compute_transform(f, x, s, **hints)`` and ``_as_integral(f, x, s)``
  59. functions. If the transform cannot be computed, raise :obj:`IntegralTransformError`.
  60. Also set ``cls._name``. For instance,
  61. >>> from sympy import LaplaceTransform
  62. >>> LaplaceTransform._name
  63. 'Laplace'
  64. Implement ``self._collapse_extra`` if your function returns more than just a
  65. number and possibly a convergence condition.
  66. """
  67. @property
  68. def function(self):
  69. """ The function to be transformed. """
  70. return self.args[0]
  71. @property
  72. def function_variable(self):
  73. """ The dependent variable of the function to be transformed. """
  74. return self.args[1]
  75. @property
  76. def transform_variable(self):
  77. """ The independent transform variable. """
  78. return self.args[2]
  79. @property
  80. def free_symbols(self):
  81. """
  82. This method returns the symbols that will exist when the transform
  83. is evaluated.
  84. """
  85. return self.function.free_symbols.union({self.transform_variable}) \
  86. - {self.function_variable}
  87. def _compute_transform(self, f, x, s, **hints):
  88. raise NotImplementedError
  89. def _as_integral(self, f, x, s):
  90. raise NotImplementedError
  91. def _collapse_extra(self, extra):
  92. cond = And(*extra)
  93. if cond == False:
  94. raise IntegralTransformError(self.__class__.name, None, '')
  95. return cond
  96. def _try_directly(self, **hints):
  97. T = None
  98. try_directly = not any(func.has(self.function_variable)
  99. for func in self.function.atoms(AppliedUndef))
  100. if try_directly:
  101. try:
  102. T = self._compute_transform(self.function,
  103. self.function_variable, self.transform_variable, **hints)
  104. except IntegralTransformError:
  105. debug('[IT _try ] Caught IntegralTransformError, returns None')
  106. T = None
  107. fn = self.function
  108. if not fn.is_Add:
  109. fn = expand_mul(fn)
  110. return fn, T
  111. def doit(self, **hints):
  112. """
  113. Try to evaluate the transform in closed form.
  114. Explanation
  115. ===========
  116. This general function handles linearity, but apart from that leaves
  117. pretty much everything to _compute_transform.
  118. Standard hints are the following:
  119. - ``simplify``: whether or not to simplify the result
  120. - ``noconds``: if True, do not return convergence conditions
  121. - ``needeval``: if True, raise IntegralTransformError instead of
  122. returning IntegralTransform objects
  123. The default values of these hints depend on the concrete transform,
  124. usually the default is
  125. ``(simplify, noconds, needeval) = (True, False, False)``.
  126. """
  127. needeval = hints.pop('needeval', False)
  128. simplify = hints.pop('simplify', True)
  129. hints['simplify'] = simplify
  130. fn, T = self._try_directly(**hints)
  131. if T is not None:
  132. return T
  133. if fn.is_Add:
  134. hints['needeval'] = needeval
  135. res = [self.__class__(*([x] + list(self.args[1:]))).doit(**hints)
  136. for x in fn.args]
  137. extra = []
  138. ress = []
  139. for x in res:
  140. if not isinstance(x, tuple):
  141. x = [x]
  142. ress.append(x[0])
  143. if len(x) == 2:
  144. # only a condition
  145. extra.append(x[1])
  146. elif len(x) > 2:
  147. # some region parameters and a condition (Mellin, Laplace)
  148. extra += [x[1:]]
  149. if simplify==True:
  150. res = Add(*ress).simplify()
  151. else:
  152. res = Add(*ress)
  153. if not extra:
  154. return res
  155. try:
  156. extra = self._collapse_extra(extra)
  157. if iterable(extra):
  158. return (res,) + tuple(extra)
  159. else:
  160. return (res, extra)
  161. except IntegralTransformError:
  162. pass
  163. if needeval:
  164. raise IntegralTransformError(
  165. self.__class__._name, self.function, 'needeval')
  166. # TODO handle derivatives etc
  167. # pull out constant coefficients
  168. coeff, rest = fn.as_coeff_mul(self.function_variable)
  169. return coeff*self.__class__(*([Mul(*rest)] + list(self.args[1:])))
  170. @property
  171. def as_integral(self):
  172. return self._as_integral(self.function, self.function_variable,
  173. self.transform_variable)
  174. def _eval_rewrite_as_Integral(self, *args, **kwargs):
  175. return self.as_integral
  176. def _simplify(expr, doit):
  177. if doit:
  178. from sympy.simplify import simplify
  179. from sympy.simplify.powsimp import powdenest
  180. return simplify(powdenest(piecewise_fold(expr), polar=True))
  181. return expr
  182. def _noconds_(default):
  183. """
  184. This is a decorator generator for dropping convergence conditions.
  185. Explanation
  186. ===========
  187. Suppose you define a function ``transform(*args)`` which returns a tuple of
  188. the form ``(result, cond1, cond2, ...)``.
  189. Decorating it ``@_noconds_(default)`` will add a new keyword argument
  190. ``noconds`` to it. If ``noconds=True``, the return value will be altered to
  191. be only ``result``, whereas if ``noconds=False`` the return value will not
  192. be altered.
  193. The default value of the ``noconds`` keyword will be ``default`` (i.e. the
  194. argument of this function).
  195. """
  196. def make_wrapper(func):
  197. @wraps(func)
  198. def wrapper(*args, noconds=default, **kwargs):
  199. res = func(*args, **kwargs)
  200. if noconds:
  201. return res[0]
  202. return res
  203. return wrapper
  204. return make_wrapper
  205. _noconds = _noconds_(False)
  206. ##########################################################################
  207. # Mellin Transform
  208. ##########################################################################
  209. def _default_integrator(f, x):
  210. return integrate(f, (x, S.Zero, S.Infinity))
  211. @_noconds
  212. def _mellin_transform(f, x, s_, integrator=_default_integrator, simplify=True):
  213. """ Backend function to compute Mellin transforms. """
  214. # We use a fresh dummy, because assumptions on s might drop conditions on
  215. # convergence of the integral.
  216. s = _dummy('s', 'mellin-transform', f)
  217. F = integrator(x**(s - 1) * f, x)
  218. if not F.has(Integral):
  219. return _simplify(F.subs(s, s_), simplify), (S.NegativeInfinity, S.Infinity), S.true
  220. if not F.is_Piecewise: # XXX can this work if integration gives continuous result now?
  221. raise IntegralTransformError('Mellin', f, 'could not compute integral')
  222. F, cond = F.args[0]
  223. if F.has(Integral):
  224. raise IntegralTransformError(
  225. 'Mellin', f, 'integral in unexpected form')
  226. def process_conds(cond):
  227. """
  228. Turn ``cond`` into a strip (a, b), and auxiliary conditions.
  229. """
  230. from sympy.solvers.inequalities import _solve_inequality
  231. a = S.NegativeInfinity
  232. b = S.Infinity
  233. aux = S.true
  234. conds = conjuncts(to_cnf(cond))
  235. t = Dummy('t', real=True)
  236. for c in conds:
  237. a_ = S.Infinity
  238. b_ = S.NegativeInfinity
  239. aux_ = []
  240. for d in disjuncts(c):
  241. d_ = d.replace(
  242. re, lambda x: x.as_real_imag()[0]).subs(re(s), t)
  243. if not d.is_Relational or \
  244. d.rel_op in ('==', '!=') \
  245. or d_.has(s) or not d_.has(t):
  246. aux_ += [d]
  247. continue
  248. soln = _solve_inequality(d_, t)
  249. if not soln.is_Relational or \
  250. soln.rel_op in ('==', '!='):
  251. aux_ += [d]
  252. continue
  253. if soln.lts == t:
  254. b_ = Max(soln.gts, b_)
  255. else:
  256. a_ = Min(soln.lts, a_)
  257. if a_ is not S.Infinity and a_ != b:
  258. a = Max(a_, a)
  259. elif b_ is not S.NegativeInfinity and b_ != a:
  260. b = Min(b_, b)
  261. else:
  262. aux = And(aux, Or(*aux_))
  263. return a, b, aux
  264. conds = [process_conds(c) for c in disjuncts(cond)]
  265. conds = [x for x in conds if x[2] != False]
  266. conds.sort(key=lambda x: (x[0] - x[1], count_ops(x[2])))
  267. if not conds:
  268. raise IntegralTransformError('Mellin', f, 'no convergence found')
  269. a, b, aux = conds[0]
  270. return _simplify(F.subs(s, s_), simplify), (a, b), aux
  271. class MellinTransform(IntegralTransform):
  272. """
  273. Class representing unevaluated Mellin transforms.
  274. For usage of this class, see the :class:`IntegralTransform` docstring.
  275. For how to compute Mellin transforms, see the :func:`mellin_transform`
  276. docstring.
  277. """
  278. _name = 'Mellin'
  279. def _compute_transform(self, f, x, s, **hints):
  280. return _mellin_transform(f, x, s, **hints)
  281. def _as_integral(self, f, x, s):
  282. return Integral(f*x**(s - 1), (x, S.Zero, S.Infinity))
  283. def _collapse_extra(self, extra):
  284. a = []
  285. b = []
  286. cond = []
  287. for (sa, sb), c in extra:
  288. a += [sa]
  289. b += [sb]
  290. cond += [c]
  291. res = (Max(*a), Min(*b)), And(*cond)
  292. if (res[0][0] >= res[0][1]) == True or res[1] == False:
  293. raise IntegralTransformError(
  294. 'Mellin', None, 'no combined convergence.')
  295. return res
  296. def mellin_transform(f, x, s, **hints):
  297. r"""
  298. Compute the Mellin transform `F(s)` of `f(x)`,
  299. .. math :: F(s) = \int_0^\infty x^{s-1} f(x) \mathrm{d}x.
  300. For all "sensible" functions, this converges absolutely in a strip
  301. `a < \operatorname{Re}(s) < b`.
  302. Explanation
  303. ===========
  304. The Mellin transform is related via change of variables to the Fourier
  305. transform, and also to the (bilateral) Laplace transform.
  306. This function returns ``(F, (a, b), cond)``
  307. where ``F`` is the Mellin transform of ``f``, ``(a, b)`` is the fundamental strip
  308. (as above), and ``cond`` are auxiliary convergence conditions.
  309. If the integral cannot be computed in closed form, this function returns
  310. an unevaluated :class:`MellinTransform` object.
  311. For a description of possible hints, refer to the docstring of
  312. :func:`sympy.integrals.transforms.IntegralTransform.doit`. If ``noconds=False``,
  313. then only `F` will be returned (i.e. not ``cond``, and also not the strip
  314. ``(a, b)``).
  315. Examples
  316. ========
  317. >>> from sympy import mellin_transform, exp
  318. >>> from sympy.abc import x, s
  319. >>> mellin_transform(exp(-x), x, s)
  320. (gamma(s), (0, oo), True)
  321. See Also
  322. ========
  323. inverse_mellin_transform, laplace_transform, fourier_transform
  324. hankel_transform, inverse_hankel_transform
  325. """
  326. return MellinTransform(f, x, s).doit(**hints)
  327. def _rewrite_sin(m_n, s, a, b):
  328. """
  329. Re-write the sine function ``sin(m*s + n)`` as gamma functions, compatible
  330. with the strip (a, b).
  331. Return ``(gamma1, gamma2, fac)`` so that ``f == fac/(gamma1 * gamma2)``.
  332. Examples
  333. ========
  334. >>> from sympy.integrals.transforms import _rewrite_sin
  335. >>> from sympy import pi, S
  336. >>> from sympy.abc import s
  337. >>> _rewrite_sin((pi, 0), s, 0, 1)
  338. (gamma(s), gamma(1 - s), pi)
  339. >>> _rewrite_sin((pi, 0), s, 1, 0)
  340. (gamma(s - 1), gamma(2 - s), -pi)
  341. >>> _rewrite_sin((pi, 0), s, -1, 0)
  342. (gamma(s + 1), gamma(-s), -pi)
  343. >>> _rewrite_sin((pi, pi/2), s, S(1)/2, S(3)/2)
  344. (gamma(s - 1/2), gamma(3/2 - s), -pi)
  345. >>> _rewrite_sin((pi, pi), s, 0, 1)
  346. (gamma(s), gamma(1 - s), -pi)
  347. >>> _rewrite_sin((2*pi, 0), s, 0, S(1)/2)
  348. (gamma(2*s), gamma(1 - 2*s), pi)
  349. >>> _rewrite_sin((2*pi, 0), s, S(1)/2, 1)
  350. (gamma(2*s - 1), gamma(2 - 2*s), -pi)
  351. """
  352. # (This is a separate function because it is moderately complicated,
  353. # and I want to doctest it.)
  354. # We want to use pi/sin(pi*x) = gamma(x)*gamma(1-x).
  355. # But there is one comlication: the gamma functions determine the
  356. # inegration contour in the definition of the G-function. Usually
  357. # it would not matter if this is slightly shifted, unless this way
  358. # we create an undefined function!
  359. # So we try to write this in such a way that the gammas are
  360. # eminently on the right side of the strip.
  361. m, n = m_n
  362. m = expand_mul(m/pi)
  363. n = expand_mul(n/pi)
  364. r = ceiling(-m*a - n.as_real_imag()[0]) # Don't use re(n), does not expand
  365. return gamma(m*s + n + r), gamma(1 - n - r - m*s), (-1)**r*pi
  366. class MellinTransformStripError(ValueError):
  367. """
  368. Exception raised by _rewrite_gamma. Mainly for internal use.
  369. """
  370. pass
  371. def _rewrite_gamma(f, s, a, b):
  372. """
  373. Try to rewrite the product f(s) as a product of gamma functions,
  374. so that the inverse Mellin transform of f can be expressed as a meijer
  375. G function.
  376. Explanation
  377. ===========
  378. Return (an, ap), (bm, bq), arg, exp, fac such that
  379. G((an, ap), (bm, bq), arg/z**exp)*fac is the inverse Mellin transform of f(s).
  380. Raises IntegralTransformError or MellinTransformStripError on failure.
  381. It is asserted that f has no poles in the fundamental strip designated by
  382. (a, b). One of a and b is allowed to be None. The fundamental strip is
  383. important, because it determines the inversion contour.
  384. This function can handle exponentials, linear factors, trigonometric
  385. functions.
  386. This is a helper function for inverse_mellin_transform that will not
  387. attempt any transformations on f.
  388. Examples
  389. ========
  390. >>> from sympy.integrals.transforms import _rewrite_gamma
  391. >>> from sympy.abc import s
  392. >>> from sympy import oo
  393. >>> _rewrite_gamma(s*(s+3)*(s-1), s, -oo, oo)
  394. (([], [-3, 0, 1]), ([-2, 1, 2], []), 1, 1, -1)
  395. >>> _rewrite_gamma((s-1)**2, s, -oo, oo)
  396. (([], [1, 1]), ([2, 2], []), 1, 1, 1)
  397. Importance of the fundamental strip:
  398. >>> _rewrite_gamma(1/s, s, 0, oo)
  399. (([1], []), ([], [0]), 1, 1, 1)
  400. >>> _rewrite_gamma(1/s, s, None, oo)
  401. (([1], []), ([], [0]), 1, 1, 1)
  402. >>> _rewrite_gamma(1/s, s, 0, None)
  403. (([1], []), ([], [0]), 1, 1, 1)
  404. >>> _rewrite_gamma(1/s, s, -oo, 0)
  405. (([], [1]), ([0], []), 1, 1, -1)
  406. >>> _rewrite_gamma(1/s, s, None, 0)
  407. (([], [1]), ([0], []), 1, 1, -1)
  408. >>> _rewrite_gamma(1/s, s, -oo, None)
  409. (([], [1]), ([0], []), 1, 1, -1)
  410. >>> _rewrite_gamma(2**(-s+3), s, -oo, oo)
  411. (([], []), ([], []), 1/2, 1, 8)
  412. """
  413. # Our strategy will be as follows:
  414. # 1) Guess a constant c such that the inversion integral should be
  415. # performed wrt s'=c*s (instead of plain s). Write s for s'.
  416. # 2) Process all factors, rewrite them independently as gamma functions in
  417. # argument s, or exponentials of s.
  418. # 3) Try to transform all gamma functions s.t. they have argument
  419. # a+s or a-s.
  420. # 4) Check that the resulting G function parameters are valid.
  421. # 5) Combine all the exponentials.
  422. a_, b_ = S([a, b])
  423. def left(c, is_numer):
  424. """
  425. Decide whether pole at c lies to the left of the fundamental strip.
  426. """
  427. # heuristically, this is the best chance for us to solve the inequalities
  428. c = expand(re(c))
  429. if a_ is None and b_ is S.Infinity:
  430. return True
  431. if a_ is None:
  432. return c < b_
  433. if b_ is None:
  434. return c <= a_
  435. if (c >= b_) == True:
  436. return False
  437. if (c <= a_) == True:
  438. return True
  439. if is_numer:
  440. return None
  441. if a_.free_symbols or b_.free_symbols or c.free_symbols:
  442. return None # XXX
  443. #raise IntegralTransformError('Inverse Mellin', f,
  444. # 'Could not determine position of singularity %s'
  445. # ' relative to fundamental strip' % c)
  446. raise MellinTransformStripError('Pole inside critical strip?')
  447. # 1)
  448. s_multipliers = []
  449. for g in f.atoms(gamma):
  450. if not g.has(s):
  451. continue
  452. arg = g.args[0]
  453. if arg.is_Add:
  454. arg = arg.as_independent(s)[1]
  455. coeff, _ = arg.as_coeff_mul(s)
  456. s_multipliers += [coeff]
  457. for g in f.atoms(sin, cos, tan, cot):
  458. if not g.has(s):
  459. continue
  460. arg = g.args[0]
  461. if arg.is_Add:
  462. arg = arg.as_independent(s)[1]
  463. coeff, _ = arg.as_coeff_mul(s)
  464. s_multipliers += [coeff/pi]
  465. s_multipliers = [Abs(x) if x.is_extended_real else x for x in s_multipliers]
  466. common_coefficient = S.One
  467. for x in s_multipliers:
  468. if not x.is_Rational:
  469. common_coefficient = x
  470. break
  471. s_multipliers = [x/common_coefficient for x in s_multipliers]
  472. if not (all(x.is_Rational for x in s_multipliers) and
  473. common_coefficient.is_extended_real):
  474. raise IntegralTransformError("Gamma", None, "Nonrational multiplier")
  475. s_multiplier = common_coefficient/reduce(ilcm, [S(x.q)
  476. for x in s_multipliers], S.One)
  477. if s_multiplier == common_coefficient:
  478. if len(s_multipliers) == 0:
  479. s_multiplier = common_coefficient
  480. else:
  481. s_multiplier = common_coefficient \
  482. *reduce(igcd, [S(x.p) for x in s_multipliers])
  483. f = f.subs(s, s/s_multiplier)
  484. fac = S.One/s_multiplier
  485. exponent = S.One/s_multiplier
  486. if a_ is not None:
  487. a_ *= s_multiplier
  488. if b_ is not None:
  489. b_ *= s_multiplier
  490. # 2)
  491. numer, denom = f.as_numer_denom()
  492. numer = Mul.make_args(numer)
  493. denom = Mul.make_args(denom)
  494. args = list(zip(numer, repeat(True))) + list(zip(denom, repeat(False)))
  495. facs = []
  496. dfacs = []
  497. # *_gammas will contain pairs (a, c) representing Gamma(a*s + c)
  498. numer_gammas = []
  499. denom_gammas = []
  500. # exponentials will contain bases for exponentials of s
  501. exponentials = []
  502. def exception(fact):
  503. return IntegralTransformError("Inverse Mellin", f, "Unrecognised form '%s'." % fact)
  504. while args:
  505. fact, is_numer = args.pop()
  506. if is_numer:
  507. ugammas, lgammas = numer_gammas, denom_gammas
  508. ufacs = facs
  509. else:
  510. ugammas, lgammas = denom_gammas, numer_gammas
  511. ufacs = dfacs
  512. def linear_arg(arg):
  513. """ Test if arg is of form a*s+b, raise exception if not. """
  514. if not arg.is_polynomial(s):
  515. raise exception(fact)
  516. p = Poly(arg, s)
  517. if p.degree() != 1:
  518. raise exception(fact)
  519. return p.all_coeffs()
  520. # constants
  521. if not fact.has(s):
  522. ufacs += [fact]
  523. # exponentials
  524. elif fact.is_Pow or isinstance(fact, exp):
  525. if fact.is_Pow:
  526. base = fact.base
  527. exp_ = fact.exp
  528. else:
  529. base = exp_polar(1)
  530. exp_ = fact.exp
  531. if exp_.is_Integer:
  532. cond = is_numer
  533. if exp_ < 0:
  534. cond = not cond
  535. args += [(base, cond)]*Abs(exp_)
  536. continue
  537. elif not base.has(s):
  538. a, b = linear_arg(exp_)
  539. if not is_numer:
  540. base = 1/base
  541. exponentials += [base**a]
  542. facs += [base**b]
  543. else:
  544. raise exception(fact)
  545. # linear factors
  546. elif fact.is_polynomial(s):
  547. p = Poly(fact, s)
  548. if p.degree() != 1:
  549. # We completely factor the poly. For this we need the roots.
  550. # Now roots() only works in some cases (low degree), and CRootOf
  551. # only works without parameters. So try both...
  552. coeff = p.LT()[1]
  553. rs = roots(p, s)
  554. if len(rs) != p.degree():
  555. rs = CRootOf.all_roots(p)
  556. ufacs += [coeff]
  557. args += [(s - c, is_numer) for c in rs]
  558. continue
  559. a, c = p.all_coeffs()
  560. ufacs += [a]
  561. c /= -a
  562. # Now need to convert s - c
  563. if left(c, is_numer):
  564. ugammas += [(S.One, -c + 1)]
  565. lgammas += [(S.One, -c)]
  566. else:
  567. ufacs += [-1]
  568. ugammas += [(S.NegativeOne, c + 1)]
  569. lgammas += [(S.NegativeOne, c)]
  570. elif isinstance(fact, gamma):
  571. a, b = linear_arg(fact.args[0])
  572. if is_numer:
  573. if (a > 0 and (left(-b/a, is_numer) == False)) or \
  574. (a < 0 and (left(-b/a, is_numer) == True)):
  575. raise NotImplementedError(
  576. 'Gammas partially over the strip.')
  577. ugammas += [(a, b)]
  578. elif isinstance(fact, sin):
  579. # We try to re-write all trigs as gammas. This is not in
  580. # general the best strategy, since sometimes this is impossible,
  581. # but rewriting as exponentials would work. However trig functions
  582. # in inverse mellin transforms usually all come from simplifying
  583. # gamma terms, so this should work.
  584. a = fact.args[0]
  585. if is_numer:
  586. # No problem with the poles.
  587. gamma1, gamma2, fac_ = gamma(a/pi), gamma(1 - a/pi), pi
  588. else:
  589. gamma1, gamma2, fac_ = _rewrite_sin(linear_arg(a), s, a_, b_)
  590. args += [(gamma1, not is_numer), (gamma2, not is_numer)]
  591. ufacs += [fac_]
  592. elif isinstance(fact, tan):
  593. a = fact.args[0]
  594. args += [(sin(a, evaluate=False), is_numer),
  595. (sin(pi/2 - a, evaluate=False), not is_numer)]
  596. elif isinstance(fact, cos):
  597. a = fact.args[0]
  598. args += [(sin(pi/2 - a, evaluate=False), is_numer)]
  599. elif isinstance(fact, cot):
  600. a = fact.args[0]
  601. args += [(sin(pi/2 - a, evaluate=False), is_numer),
  602. (sin(a, evaluate=False), not is_numer)]
  603. else:
  604. raise exception(fact)
  605. fac *= Mul(*facs)/Mul(*dfacs)
  606. # 3)
  607. an, ap, bm, bq = [], [], [], []
  608. for gammas, plus, minus, is_numer in [(numer_gammas, an, bm, True),
  609. (denom_gammas, bq, ap, False)]:
  610. while gammas:
  611. a, c = gammas.pop()
  612. if a != -1 and a != +1:
  613. # We use the gamma function multiplication theorem.
  614. p = Abs(S(a))
  615. newa = a/p
  616. newc = c/p
  617. if not a.is_Integer:
  618. raise TypeError("a is not an integer")
  619. for k in range(p):
  620. gammas += [(newa, newc + k/p)]
  621. if is_numer:
  622. fac *= (2*pi)**((1 - p)/2) * p**(c - S.Half)
  623. exponentials += [p**a]
  624. else:
  625. fac /= (2*pi)**((1 - p)/2) * p**(c - S.Half)
  626. exponentials += [p**(-a)]
  627. continue
  628. if a == +1:
  629. plus.append(1 - c)
  630. else:
  631. minus.append(c)
  632. # 4)
  633. # TODO
  634. # 5)
  635. arg = Mul(*exponentials)
  636. # for testability, sort the arguments
  637. an.sort(key=default_sort_key)
  638. ap.sort(key=default_sort_key)
  639. bm.sort(key=default_sort_key)
  640. bq.sort(key=default_sort_key)
  641. return (an, ap), (bm, bq), arg, exponent, fac
  642. @_noconds_(True)
  643. def _inverse_mellin_transform(F, s, x_, strip, as_meijerg=False):
  644. """ A helper for the real inverse_mellin_transform function, this one here
  645. assumes x to be real and positive. """
  646. x = _dummy('t', 'inverse-mellin-transform', F, positive=True)
  647. # Actually, we won't try integration at all. Instead we use the definition
  648. # of the Meijer G function as a fairly general inverse mellin transform.
  649. F = F.rewrite(gamma)
  650. for g in [factor(F), expand_mul(F), expand(F)]:
  651. if g.is_Add:
  652. # do all terms separately
  653. ress = [_inverse_mellin_transform(G, s, x, strip, as_meijerg,
  654. noconds=False)
  655. for G in g.args]
  656. conds = [p[1] for p in ress]
  657. ress = [p[0] for p in ress]
  658. res = Add(*ress)
  659. if not as_meijerg:
  660. res = factor(res, gens=res.atoms(Heaviside))
  661. return res.subs(x, x_), And(*conds)
  662. try:
  663. a, b, C, e, fac = _rewrite_gamma(g, s, strip[0], strip[1])
  664. except IntegralTransformError:
  665. continue
  666. try:
  667. G = meijerg(a, b, C/x**e)
  668. except ValueError:
  669. continue
  670. if as_meijerg:
  671. h = G
  672. else:
  673. try:
  674. from sympy.simplify import hyperexpand
  675. h = hyperexpand(G)
  676. except NotImplementedError:
  677. raise IntegralTransformError(
  678. 'Inverse Mellin', F, 'Could not calculate integral')
  679. if h.is_Piecewise and len(h.args) == 3:
  680. # XXX we break modularity here!
  681. h = Heaviside(x - Abs(C))*h.args[0].args[0] \
  682. + Heaviside(Abs(C) - x)*h.args[1].args[0]
  683. # We must ensure that the integral along the line we want converges,
  684. # and return that value.
  685. # See [L], 5.2
  686. cond = [Abs(arg(G.argument)) < G.delta*pi]
  687. # Note: we allow ">=" here, this corresponds to convergence if we let
  688. # limits go to oo symmetrically. ">" corresponds to absolute convergence.
  689. cond += [And(Or(len(G.ap) != len(G.bq), 0 >= re(G.nu) + 1),
  690. Abs(arg(G.argument)) == G.delta*pi)]
  691. cond = Or(*cond)
  692. if cond == False:
  693. raise IntegralTransformError(
  694. 'Inverse Mellin', F, 'does not converge')
  695. return (h*fac).subs(x, x_), cond
  696. raise IntegralTransformError('Inverse Mellin', F, '')
  697. _allowed = None
  698. class InverseMellinTransform(IntegralTransform):
  699. """
  700. Class representing unevaluated inverse Mellin transforms.
  701. For usage of this class, see the :class:`IntegralTransform` docstring.
  702. For how to compute inverse Mellin transforms, see the
  703. :func:`inverse_mellin_transform` docstring.
  704. """
  705. _name = 'Inverse Mellin'
  706. _none_sentinel = Dummy('None')
  707. _c = Dummy('c')
  708. def __new__(cls, F, s, x, a, b, **opts):
  709. if a is None:
  710. a = InverseMellinTransform._none_sentinel
  711. if b is None:
  712. b = InverseMellinTransform._none_sentinel
  713. return IntegralTransform.__new__(cls, F, s, x, a, b, **opts)
  714. @property
  715. def fundamental_strip(self):
  716. a, b = self.args[3], self.args[4]
  717. if a is InverseMellinTransform._none_sentinel:
  718. a = None
  719. if b is InverseMellinTransform._none_sentinel:
  720. b = None
  721. return a, b
  722. def _compute_transform(self, F, s, x, **hints):
  723. # IntegralTransform's doit will cause this hint to exist, but
  724. # InverseMellinTransform should ignore it
  725. hints.pop('simplify', True)
  726. global _allowed
  727. if _allowed is None:
  728. _allowed = {
  729. exp, gamma, sin, cos, tan, cot, cosh, sinh, tanh, coth,
  730. factorial, rf}
  731. for f in postorder_traversal(F):
  732. if f.is_Function and f.has(s) and f.func not in _allowed:
  733. raise IntegralTransformError('Inverse Mellin', F,
  734. 'Component %s not recognised.' % f)
  735. strip = self.fundamental_strip
  736. return _inverse_mellin_transform(F, s, x, strip, **hints)
  737. def _as_integral(self, F, s, x):
  738. c = self.__class__._c
  739. return Integral(F*x**(-s), (s, c - S.ImaginaryUnit*S.Infinity, c +
  740. S.ImaginaryUnit*S.Infinity))/(2*S.Pi*S.ImaginaryUnit)
  741. def inverse_mellin_transform(F, s, x, strip, **hints):
  742. r"""
  743. Compute the inverse Mellin transform of `F(s)` over the fundamental
  744. strip given by ``strip=(a, b)``.
  745. Explanation
  746. ===========
  747. This can be defined as
  748. .. math:: f(x) = \frac{1}{2\pi i} \int_{c - i\infty}^{c + i\infty} x^{-s} F(s) \mathrm{d}s,
  749. for any `c` in the fundamental strip. Under certain regularity
  750. conditions on `F` and/or `f`,
  751. this recovers `f` from its Mellin transform `F`
  752. (and vice versa), for positive real `x`.
  753. One of `a` or `b` may be passed as ``None``; a suitable `c` will be
  754. inferred.
  755. If the integral cannot be computed in closed form, this function returns
  756. an unevaluated :class:`InverseMellinTransform` object.
  757. Note that this function will assume x to be positive and real, regardless
  758. of the SymPy assumptions!
  759. For a description of possible hints, refer to the docstring of
  760. :func:`sympy.integrals.transforms.IntegralTransform.doit`.
  761. Examples
  762. ========
  763. >>> from sympy import inverse_mellin_transform, oo, gamma
  764. >>> from sympy.abc import x, s
  765. >>> inverse_mellin_transform(gamma(s), s, x, (0, oo))
  766. exp(-x)
  767. The fundamental strip matters:
  768. >>> f = 1/(s**2 - 1)
  769. >>> inverse_mellin_transform(f, s, x, (-oo, -1))
  770. x*(1 - 1/x**2)*Heaviside(x - 1)/2
  771. >>> inverse_mellin_transform(f, s, x, (-1, 1))
  772. -x*Heaviside(1 - x)/2 - Heaviside(x - 1)/(2*x)
  773. >>> inverse_mellin_transform(f, s, x, (1, oo))
  774. (1/2 - x**2/2)*Heaviside(1 - x)/x
  775. See Also
  776. ========
  777. mellin_transform
  778. hankel_transform, inverse_hankel_transform
  779. """
  780. return InverseMellinTransform(F, s, x, strip[0], strip[1]).doit(**hints)
  781. ##########################################################################
  782. # Fourier Transform
  783. ##########################################################################
  784. @_noconds_(True)
  785. def _fourier_transform(f, x, k, a, b, name, simplify=True):
  786. r"""
  787. Compute a general Fourier-type transform
  788. .. math::
  789. F(k) = a \int_{-\infty}^{\infty} e^{bixk} f(x)\, dx.
  790. For suitable choice of *a* and *b*, this reduces to the standard Fourier
  791. and inverse Fourier transforms.
  792. """
  793. F = integrate(a*f*exp(b*S.ImaginaryUnit*x*k), (x, S.NegativeInfinity, S.Infinity))
  794. if not F.has(Integral):
  795. return _simplify(F, simplify), S.true
  796. integral_f = integrate(f, (x, S.NegativeInfinity, S.Infinity))
  797. if integral_f in (S.NegativeInfinity, S.Infinity, S.NaN) or integral_f.has(Integral):
  798. raise IntegralTransformError(name, f, 'function not integrable on real axis')
  799. if not F.is_Piecewise:
  800. raise IntegralTransformError(name, f, 'could not compute integral')
  801. F, cond = F.args[0]
  802. if F.has(Integral):
  803. raise IntegralTransformError(name, f, 'integral in unexpected form')
  804. return _simplify(F, simplify), cond
  805. class FourierTypeTransform(IntegralTransform):
  806. """ Base class for Fourier transforms."""
  807. def a(self):
  808. raise NotImplementedError(
  809. "Class %s must implement a(self) but does not" % self.__class__)
  810. def b(self):
  811. raise NotImplementedError(
  812. "Class %s must implement b(self) but does not" % self.__class__)
  813. def _compute_transform(self, f, x, k, **hints):
  814. return _fourier_transform(f, x, k,
  815. self.a(), self.b(),
  816. self.__class__._name, **hints)
  817. def _as_integral(self, f, x, k):
  818. a = self.a()
  819. b = self.b()
  820. return Integral(a*f*exp(b*S.ImaginaryUnit*x*k), (x, S.NegativeInfinity, S.Infinity))
  821. class FourierTransform(FourierTypeTransform):
  822. """
  823. Class representing unevaluated Fourier transforms.
  824. For usage of this class, see the :class:`IntegralTransform` docstring.
  825. For how to compute Fourier transforms, see the :func:`fourier_transform`
  826. docstring.
  827. """
  828. _name = 'Fourier'
  829. def a(self):
  830. return 1
  831. def b(self):
  832. return -2*S.Pi
  833. def fourier_transform(f, x, k, **hints):
  834. r"""
  835. Compute the unitary, ordinary-frequency Fourier transform of ``f``, defined
  836. as
  837. .. math:: F(k) = \int_{-\infty}^\infty f(x) e^{-2\pi i x k} \mathrm{d} x.
  838. Explanation
  839. ===========
  840. If the transform cannot be computed in closed form, this
  841. function returns an unevaluated :class:`FourierTransform` object.
  842. For other Fourier transform conventions, see the function
  843. :func:`sympy.integrals.transforms._fourier_transform`.
  844. For a description of possible hints, refer to the docstring of
  845. :func:`sympy.integrals.transforms.IntegralTransform.doit`.
  846. Note that for this transform, by default ``noconds=True``.
  847. Examples
  848. ========
  849. >>> from sympy import fourier_transform, exp
  850. >>> from sympy.abc import x, k
  851. >>> fourier_transform(exp(-x**2), x, k)
  852. sqrt(pi)*exp(-pi**2*k**2)
  853. >>> fourier_transform(exp(-x**2), x, k, noconds=False)
  854. (sqrt(pi)*exp(-pi**2*k**2), True)
  855. See Also
  856. ========
  857. inverse_fourier_transform
  858. sine_transform, inverse_sine_transform
  859. cosine_transform, inverse_cosine_transform
  860. hankel_transform, inverse_hankel_transform
  861. mellin_transform, laplace_transform
  862. """
  863. return FourierTransform(f, x, k).doit(**hints)
  864. class InverseFourierTransform(FourierTypeTransform):
  865. """
  866. Class representing unevaluated inverse Fourier transforms.
  867. For usage of this class, see the :class:`IntegralTransform` docstring.
  868. For how to compute inverse Fourier transforms, see the
  869. :func:`inverse_fourier_transform` docstring.
  870. """
  871. _name = 'Inverse Fourier'
  872. def a(self):
  873. return 1
  874. def b(self):
  875. return 2*S.Pi
  876. def inverse_fourier_transform(F, k, x, **hints):
  877. r"""
  878. Compute the unitary, ordinary-frequency inverse Fourier transform of `F`,
  879. defined as
  880. .. math:: f(x) = \int_{-\infty}^\infty F(k) e^{2\pi i x k} \mathrm{d} k.
  881. Explanation
  882. ===========
  883. If the transform cannot be computed in closed form, this
  884. function returns an unevaluated :class:`InverseFourierTransform` object.
  885. For other Fourier transform conventions, see the function
  886. :func:`sympy.integrals.transforms._fourier_transform`.
  887. For a description of possible hints, refer to the docstring of
  888. :func:`sympy.integrals.transforms.IntegralTransform.doit`.
  889. Note that for this transform, by default ``noconds=True``.
  890. Examples
  891. ========
  892. >>> from sympy import inverse_fourier_transform, exp, sqrt, pi
  893. >>> from sympy.abc import x, k
  894. >>> inverse_fourier_transform(sqrt(pi)*exp(-(pi*k)**2), k, x)
  895. exp(-x**2)
  896. >>> inverse_fourier_transform(sqrt(pi)*exp(-(pi*k)**2), k, x, noconds=False)
  897. (exp(-x**2), True)
  898. See Also
  899. ========
  900. fourier_transform
  901. sine_transform, inverse_sine_transform
  902. cosine_transform, inverse_cosine_transform
  903. hankel_transform, inverse_hankel_transform
  904. mellin_transform, laplace_transform
  905. """
  906. return InverseFourierTransform(F, k, x).doit(**hints)
  907. ##########################################################################
  908. # Fourier Sine and Cosine Transform
  909. ##########################################################################
  910. @_noconds_(True)
  911. def _sine_cosine_transform(f, x, k, a, b, K, name, simplify=True):
  912. """
  913. Compute a general sine or cosine-type transform
  914. F(k) = a int_0^oo b*sin(x*k) f(x) dx.
  915. F(k) = a int_0^oo b*cos(x*k) f(x) dx.
  916. For suitable choice of a and b, this reduces to the standard sine/cosine
  917. and inverse sine/cosine transforms.
  918. """
  919. F = integrate(a*f*K(b*x*k), (x, S.Zero, S.Infinity))
  920. if not F.has(Integral):
  921. return _simplify(F, simplify), S.true
  922. if not F.is_Piecewise:
  923. raise IntegralTransformError(name, f, 'could not compute integral')
  924. F, cond = F.args[0]
  925. if F.has(Integral):
  926. raise IntegralTransformError(name, f, 'integral in unexpected form')
  927. return _simplify(F, simplify), cond
  928. class SineCosineTypeTransform(IntegralTransform):
  929. """
  930. Base class for sine and cosine transforms.
  931. Specify cls._kern.
  932. """
  933. def a(self):
  934. raise NotImplementedError(
  935. "Class %s must implement a(self) but does not" % self.__class__)
  936. def b(self):
  937. raise NotImplementedError(
  938. "Class %s must implement b(self) but does not" % self.__class__)
  939. def _compute_transform(self, f, x, k, **hints):
  940. return _sine_cosine_transform(f, x, k,
  941. self.a(), self.b(),
  942. self.__class__._kern,
  943. self.__class__._name, **hints)
  944. def _as_integral(self, f, x, k):
  945. a = self.a()
  946. b = self.b()
  947. K = self.__class__._kern
  948. return Integral(a*f*K(b*x*k), (x, S.Zero, S.Infinity))
  949. class SineTransform(SineCosineTypeTransform):
  950. """
  951. Class representing unevaluated sine transforms.
  952. For usage of this class, see the :class:`IntegralTransform` docstring.
  953. For how to compute sine transforms, see the :func:`sine_transform`
  954. docstring.
  955. """
  956. _name = 'Sine'
  957. _kern = sin
  958. def a(self):
  959. return sqrt(2)/sqrt(pi)
  960. def b(self):
  961. return S.One
  962. def sine_transform(f, x, k, **hints):
  963. r"""
  964. Compute the unitary, ordinary-frequency sine transform of `f`, defined
  965. as
  966. .. math:: F(k) = \sqrt{\frac{2}{\pi}} \int_{0}^\infty f(x) \sin(2\pi x k) \mathrm{d} x.
  967. Explanation
  968. ===========
  969. If the transform cannot be computed in closed form, this
  970. function returns an unevaluated :class:`SineTransform` object.
  971. For a description of possible hints, refer to the docstring of
  972. :func:`sympy.integrals.transforms.IntegralTransform.doit`.
  973. Note that for this transform, by default ``noconds=True``.
  974. Examples
  975. ========
  976. >>> from sympy import sine_transform, exp
  977. >>> from sympy.abc import x, k, a
  978. >>> sine_transform(x*exp(-a*x**2), x, k)
  979. sqrt(2)*k*exp(-k**2/(4*a))/(4*a**(3/2))
  980. >>> sine_transform(x**(-a), x, k)
  981. 2**(1/2 - a)*k**(a - 1)*gamma(1 - a/2)/gamma(a/2 + 1/2)
  982. See Also
  983. ========
  984. fourier_transform, inverse_fourier_transform
  985. inverse_sine_transform
  986. cosine_transform, inverse_cosine_transform
  987. hankel_transform, inverse_hankel_transform
  988. mellin_transform, laplace_transform
  989. """
  990. return SineTransform(f, x, k).doit(**hints)
  991. class InverseSineTransform(SineCosineTypeTransform):
  992. """
  993. Class representing unevaluated inverse sine transforms.
  994. For usage of this class, see the :class:`IntegralTransform` docstring.
  995. For how to compute inverse sine transforms, see the
  996. :func:`inverse_sine_transform` docstring.
  997. """
  998. _name = 'Inverse Sine'
  999. _kern = sin
  1000. def a(self):
  1001. return sqrt(2)/sqrt(pi)
  1002. def b(self):
  1003. return S.One
  1004. def inverse_sine_transform(F, k, x, **hints):
  1005. r"""
  1006. Compute the unitary, ordinary-frequency inverse sine transform of `F`,
  1007. defined as
  1008. .. math:: f(x) = \sqrt{\frac{2}{\pi}} \int_{0}^\infty F(k) \sin(2\pi x k) \mathrm{d} k.
  1009. Explanation
  1010. ===========
  1011. If the transform cannot be computed in closed form, this
  1012. function returns an unevaluated :class:`InverseSineTransform` object.
  1013. For a description of possible hints, refer to the docstring of
  1014. :func:`sympy.integrals.transforms.IntegralTransform.doit`.
  1015. Note that for this transform, by default ``noconds=True``.
  1016. Examples
  1017. ========
  1018. >>> from sympy import inverse_sine_transform, exp, sqrt, gamma
  1019. >>> from sympy.abc import x, k, a
  1020. >>> inverse_sine_transform(2**((1-2*a)/2)*k**(a - 1)*
  1021. ... gamma(-a/2 + 1)/gamma((a+1)/2), k, x)
  1022. x**(-a)
  1023. >>> inverse_sine_transform(sqrt(2)*k*exp(-k**2/(4*a))/(4*sqrt(a)**3), k, x)
  1024. x*exp(-a*x**2)
  1025. See Also
  1026. ========
  1027. fourier_transform, inverse_fourier_transform
  1028. sine_transform
  1029. cosine_transform, inverse_cosine_transform
  1030. hankel_transform, inverse_hankel_transform
  1031. mellin_transform, laplace_transform
  1032. """
  1033. return InverseSineTransform(F, k, x).doit(**hints)
  1034. class CosineTransform(SineCosineTypeTransform):
  1035. """
  1036. Class representing unevaluated cosine transforms.
  1037. For usage of this class, see the :class:`IntegralTransform` docstring.
  1038. For how to compute cosine transforms, see the :func:`cosine_transform`
  1039. docstring.
  1040. """
  1041. _name = 'Cosine'
  1042. _kern = cos
  1043. def a(self):
  1044. return sqrt(2)/sqrt(pi)
  1045. def b(self):
  1046. return S.One
  1047. def cosine_transform(f, x, k, **hints):
  1048. r"""
  1049. Compute the unitary, ordinary-frequency cosine transform of `f`, defined
  1050. as
  1051. .. math:: F(k) = \sqrt{\frac{2}{\pi}} \int_{0}^\infty f(x) \cos(2\pi x k) \mathrm{d} x.
  1052. Explanation
  1053. ===========
  1054. If the transform cannot be computed in closed form, this
  1055. function returns an unevaluated :class:`CosineTransform` object.
  1056. For a description of possible hints, refer to the docstring of
  1057. :func:`sympy.integrals.transforms.IntegralTransform.doit`.
  1058. Note that for this transform, by default ``noconds=True``.
  1059. Examples
  1060. ========
  1061. >>> from sympy import cosine_transform, exp, sqrt, cos
  1062. >>> from sympy.abc import x, k, a
  1063. >>> cosine_transform(exp(-a*x), x, k)
  1064. sqrt(2)*a/(sqrt(pi)*(a**2 + k**2))
  1065. >>> cosine_transform(exp(-a*sqrt(x))*cos(a*sqrt(x)), x, k)
  1066. a*exp(-a**2/(2*k))/(2*k**(3/2))
  1067. See Also
  1068. ========
  1069. fourier_transform, inverse_fourier_transform,
  1070. sine_transform, inverse_sine_transform
  1071. inverse_cosine_transform
  1072. hankel_transform, inverse_hankel_transform
  1073. mellin_transform, laplace_transform
  1074. """
  1075. return CosineTransform(f, x, k).doit(**hints)
  1076. class InverseCosineTransform(SineCosineTypeTransform):
  1077. """
  1078. Class representing unevaluated inverse cosine transforms.
  1079. For usage of this class, see the :class:`IntegralTransform` docstring.
  1080. For how to compute inverse cosine transforms, see the
  1081. :func:`inverse_cosine_transform` docstring.
  1082. """
  1083. _name = 'Inverse Cosine'
  1084. _kern = cos
  1085. def a(self):
  1086. return sqrt(2)/sqrt(pi)
  1087. def b(self):
  1088. return S.One
  1089. def inverse_cosine_transform(F, k, x, **hints):
  1090. r"""
  1091. Compute the unitary, ordinary-frequency inverse cosine transform of `F`,
  1092. defined as
  1093. .. math:: f(x) = \sqrt{\frac{2}{\pi}} \int_{0}^\infty F(k) \cos(2\pi x k) \mathrm{d} k.
  1094. Explanation
  1095. ===========
  1096. If the transform cannot be computed in closed form, this
  1097. function returns an unevaluated :class:`InverseCosineTransform` object.
  1098. For a description of possible hints, refer to the docstring of
  1099. :func:`sympy.integrals.transforms.IntegralTransform.doit`.
  1100. Note that for this transform, by default ``noconds=True``.
  1101. Examples
  1102. ========
  1103. >>> from sympy import inverse_cosine_transform, sqrt, pi
  1104. >>> from sympy.abc import x, k, a
  1105. >>> inverse_cosine_transform(sqrt(2)*a/(sqrt(pi)*(a**2 + k**2)), k, x)
  1106. exp(-a*x)
  1107. >>> inverse_cosine_transform(1/sqrt(k), k, x)
  1108. 1/sqrt(x)
  1109. See Also
  1110. ========
  1111. fourier_transform, inverse_fourier_transform,
  1112. sine_transform, inverse_sine_transform
  1113. cosine_transform
  1114. hankel_transform, inverse_hankel_transform
  1115. mellin_transform, laplace_transform
  1116. """
  1117. return InverseCosineTransform(F, k, x).doit(**hints)
  1118. ##########################################################################
  1119. # Hankel Transform
  1120. ##########################################################################
  1121. @_noconds_(True)
  1122. def _hankel_transform(f, r, k, nu, name, simplify=True):
  1123. r"""
  1124. Compute a general Hankel transform
  1125. .. math:: F_\nu(k) = \int_{0}^\infty f(r) J_\nu(k r) r \mathrm{d} r.
  1126. """
  1127. F = integrate(f*besselj(nu, k*r)*r, (r, S.Zero, S.Infinity))
  1128. if not F.has(Integral):
  1129. return _simplify(F, simplify), S.true
  1130. if not F.is_Piecewise:
  1131. raise IntegralTransformError(name, f, 'could not compute integral')
  1132. F, cond = F.args[0]
  1133. if F.has(Integral):
  1134. raise IntegralTransformError(name, f, 'integral in unexpected form')
  1135. return _simplify(F, simplify), cond
  1136. class HankelTypeTransform(IntegralTransform):
  1137. """
  1138. Base class for Hankel transforms.
  1139. """
  1140. def doit(self, **hints):
  1141. return self._compute_transform(self.function,
  1142. self.function_variable,
  1143. self.transform_variable,
  1144. self.args[3],
  1145. **hints)
  1146. def _compute_transform(self, f, r, k, nu, **hints):
  1147. return _hankel_transform(f, r, k, nu, self._name, **hints)
  1148. def _as_integral(self, f, r, k, nu):
  1149. return Integral(f*besselj(nu, k*r)*r, (r, S.Zero, S.Infinity))
  1150. @property
  1151. def as_integral(self):
  1152. return self._as_integral(self.function,
  1153. self.function_variable,
  1154. self.transform_variable,
  1155. self.args[3])
  1156. class HankelTransform(HankelTypeTransform):
  1157. """
  1158. Class representing unevaluated Hankel transforms.
  1159. For usage of this class, see the :class:`IntegralTransform` docstring.
  1160. For how to compute Hankel transforms, see the :func:`hankel_transform`
  1161. docstring.
  1162. """
  1163. _name = 'Hankel'
  1164. def hankel_transform(f, r, k, nu, **hints):
  1165. r"""
  1166. Compute the Hankel transform of `f`, defined as
  1167. .. math:: F_\nu(k) = \int_{0}^\infty f(r) J_\nu(k r) r \mathrm{d} r.
  1168. Explanation
  1169. ===========
  1170. If the transform cannot be computed in closed form, this
  1171. function returns an unevaluated :class:`HankelTransform` object.
  1172. For a description of possible hints, refer to the docstring of
  1173. :func:`sympy.integrals.transforms.IntegralTransform.doit`.
  1174. Note that for this transform, by default ``noconds=True``.
  1175. Examples
  1176. ========
  1177. >>> from sympy import hankel_transform, inverse_hankel_transform
  1178. >>> from sympy import exp
  1179. >>> from sympy.abc import r, k, m, nu, a
  1180. >>> ht = hankel_transform(1/r**m, r, k, nu)
  1181. >>> ht
  1182. 2*k**(m - 2)*gamma(-m/2 + nu/2 + 1)/(2**m*gamma(m/2 + nu/2))
  1183. >>> inverse_hankel_transform(ht, k, r, nu)
  1184. r**(-m)
  1185. >>> ht = hankel_transform(exp(-a*r), r, k, 0)
  1186. >>> ht
  1187. a/(k**3*(a**2/k**2 + 1)**(3/2))
  1188. >>> inverse_hankel_transform(ht, k, r, 0)
  1189. exp(-a*r)
  1190. See Also
  1191. ========
  1192. fourier_transform, inverse_fourier_transform
  1193. sine_transform, inverse_sine_transform
  1194. cosine_transform, inverse_cosine_transform
  1195. inverse_hankel_transform
  1196. mellin_transform, laplace_transform
  1197. """
  1198. return HankelTransform(f, r, k, nu).doit(**hints)
  1199. class InverseHankelTransform(HankelTypeTransform):
  1200. """
  1201. Class representing unevaluated inverse Hankel transforms.
  1202. For usage of this class, see the :class:`IntegralTransform` docstring.
  1203. For how to compute inverse Hankel transforms, see the
  1204. :func:`inverse_hankel_transform` docstring.
  1205. """
  1206. _name = 'Inverse Hankel'
  1207. def inverse_hankel_transform(F, k, r, nu, **hints):
  1208. r"""
  1209. Compute the inverse Hankel transform of `F` defined as
  1210. .. math:: f(r) = \int_{0}^\infty F_\nu(k) J_\nu(k r) k \mathrm{d} k.
  1211. Explanation
  1212. ===========
  1213. If the transform cannot be computed in closed form, this
  1214. function returns an unevaluated :class:`InverseHankelTransform` object.
  1215. For a description of possible hints, refer to the docstring of
  1216. :func:`sympy.integrals.transforms.IntegralTransform.doit`.
  1217. Note that for this transform, by default ``noconds=True``.
  1218. Examples
  1219. ========
  1220. >>> from sympy import hankel_transform, inverse_hankel_transform
  1221. >>> from sympy import exp
  1222. >>> from sympy.abc import r, k, m, nu, a
  1223. >>> ht = hankel_transform(1/r**m, r, k, nu)
  1224. >>> ht
  1225. 2*k**(m - 2)*gamma(-m/2 + nu/2 + 1)/(2**m*gamma(m/2 + nu/2))
  1226. >>> inverse_hankel_transform(ht, k, r, nu)
  1227. r**(-m)
  1228. >>> ht = hankel_transform(exp(-a*r), r, k, 0)
  1229. >>> ht
  1230. a/(k**3*(a**2/k**2 + 1)**(3/2))
  1231. >>> inverse_hankel_transform(ht, k, r, 0)
  1232. exp(-a*r)
  1233. See Also
  1234. ========
  1235. fourier_transform, inverse_fourier_transform
  1236. sine_transform, inverse_sine_transform
  1237. cosine_transform, inverse_cosine_transform
  1238. hankel_transform
  1239. mellin_transform, laplace_transform
  1240. """
  1241. return InverseHankelTransform(F, k, r, nu).doit(**hints)
  1242. ##########################################################################
  1243. # Laplace Transform
  1244. ##########################################################################
  1245. # Stub classes and functions that used to be here
  1246. import sympy.integrals.laplace as _laplace
  1247. LaplaceTransform = _laplace.LaplaceTransform
  1248. laplace_transform = _laplace.laplace_transform
  1249. InverseLaplaceTransform = _laplace.InverseLaplaceTransform
  1250. inverse_laplace_transform = _laplace.inverse_laplace_transform