Coverage for /usr/lib/python3/dist-packages/sympy/integrals/transforms.py: 26%
581 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1""" Integral Transforms """
2from functools import reduce, wraps
3from itertools import repeat
4from sympy.core import S, pi
5from sympy.core.add import Add
6from sympy.core.function import (
7 AppliedUndef, count_ops, expand, expand_mul, Function)
8from sympy.core.mul import Mul
9from sympy.core.numbers import igcd, ilcm
10from sympy.core.sorting import default_sort_key
11from sympy.core.symbol import Dummy
12from sympy.core.traversal import postorder_traversal
13from sympy.functions.combinatorial.factorials import factorial, rf
14from sympy.functions.elementary.complexes import re, arg, Abs
15from sympy.functions.elementary.exponential import exp, exp_polar
16from sympy.functions.elementary.hyperbolic import cosh, coth, sinh, tanh
17from sympy.functions.elementary.integers import ceiling
18from sympy.functions.elementary.miscellaneous import Max, Min, sqrt
19from sympy.functions.elementary.piecewise import piecewise_fold
20from sympy.functions.elementary.trigonometric import cos, cot, sin, tan
21from sympy.functions.special.bessel import besselj
22from sympy.functions.special.delta_functions import Heaviside
23from sympy.functions.special.gamma_functions import gamma
24from sympy.functions.special.hyper import meijerg
25from sympy.integrals import integrate, Integral
26from sympy.integrals.meijerint import _dummy
27from sympy.logic.boolalg import to_cnf, conjuncts, disjuncts, Or, And
28from sympy.polys.polyroots import roots
29from sympy.polys.polytools import factor, Poly
30from sympy.polys.rootoftools import CRootOf
31from sympy.utilities.iterables import iterable
32from sympy.utilities.misc import debug
35##########################################################################
36# Helpers / Utilities
37##########################################################################
40class IntegralTransformError(NotImplementedError):
41 """
42 Exception raised in relation to problems computing transforms.
44 Explanation
45 ===========
47 This class is mostly used internally; if integrals cannot be computed
48 objects representing unevaluated transforms are usually returned.
50 The hint ``needeval=True`` can be used to disable returning transform
51 objects, and instead raise this exception if an integral cannot be
52 computed.
53 """
54 def __init__(self, transform, function, msg):
55 super().__init__(
56 "%s Transform could not be computed: %s." % (transform, msg))
57 self.function = function
60class IntegralTransform(Function):
61 """
62 Base class for integral transforms.
64 Explanation
65 ===========
67 This class represents unevaluated transforms.
69 To implement a concrete transform, derive from this class and implement
70 the ``_compute_transform(f, x, s, **hints)`` and ``_as_integral(f, x, s)``
71 functions. If the transform cannot be computed, raise :obj:`IntegralTransformError`.
73 Also set ``cls._name``. For instance,
75 >>> from sympy import LaplaceTransform
76 >>> LaplaceTransform._name
77 'Laplace'
79 Implement ``self._collapse_extra`` if your function returns more than just a
80 number and possibly a convergence condition.
81 """
83 @property
84 def function(self):
85 """ The function to be transformed. """
86 return self.args[0]
88 @property
89 def function_variable(self):
90 """ The dependent variable of the function to be transformed. """
91 return self.args[1]
93 @property
94 def transform_variable(self):
95 """ The independent transform variable. """
96 return self.args[2]
98 @property
99 def free_symbols(self):
100 """
101 This method returns the symbols that will exist when the transform
102 is evaluated.
103 """
104 return self.function.free_symbols.union({self.transform_variable}) \
105 - {self.function_variable}
107 def _compute_transform(self, f, x, s, **hints):
108 raise NotImplementedError
110 def _as_integral(self, f, x, s):
111 raise NotImplementedError
113 def _collapse_extra(self, extra):
114 cond = And(*extra)
115 if cond == False:
116 raise IntegralTransformError(self.__class__.name, None, '')
117 return cond
119 def _try_directly(self, **hints):
120 T = None
121 try_directly = not any(func.has(self.function_variable)
122 for func in self.function.atoms(AppliedUndef))
123 if try_directly:
124 try:
125 T = self._compute_transform(self.function,
126 self.function_variable, self.transform_variable, **hints)
127 except IntegralTransformError:
128 debug('[IT _try ] Caught IntegralTransformError, returns None')
129 T = None
131 fn = self.function
132 if not fn.is_Add:
133 fn = expand_mul(fn)
134 return fn, T
136 def doit(self, **hints):
137 """
138 Try to evaluate the transform in closed form.
140 Explanation
141 ===========
143 This general function handles linearity, but apart from that leaves
144 pretty much everything to _compute_transform.
146 Standard hints are the following:
148 - ``simplify``: whether or not to simplify the result
149 - ``noconds``: if True, do not return convergence conditions
150 - ``needeval``: if True, raise IntegralTransformError instead of
151 returning IntegralTransform objects
153 The default values of these hints depend on the concrete transform,
154 usually the default is
155 ``(simplify, noconds, needeval) = (True, False, False)``.
156 """
157 needeval = hints.pop('needeval', False)
158 simplify = hints.pop('simplify', True)
159 hints['simplify'] = simplify
161 fn, T = self._try_directly(**hints)
163 if T is not None:
164 return T
166 if fn.is_Add:
167 hints['needeval'] = needeval
168 res = [self.__class__(*([x] + list(self.args[1:]))).doit(**hints)
169 for x in fn.args]
170 extra = []
171 ress = []
172 for x in res:
173 if not isinstance(x, tuple):
174 x = [x]
175 ress.append(x[0])
176 if len(x) == 2:
177 # only a condition
178 extra.append(x[1])
179 elif len(x) > 2:
180 # some region parameters and a condition (Mellin, Laplace)
181 extra += [x[1:]]
182 if simplify==True:
183 res = Add(*ress).simplify()
184 else:
185 res = Add(*ress)
186 if not extra:
187 return res
188 try:
189 extra = self._collapse_extra(extra)
190 if iterable(extra):
191 return (res,) + tuple(extra)
192 else:
193 return (res, extra)
194 except IntegralTransformError:
195 pass
197 if needeval:
198 raise IntegralTransformError(
199 self.__class__._name, self.function, 'needeval')
201 # TODO handle derivatives etc
203 # pull out constant coefficients
204 coeff, rest = fn.as_coeff_mul(self.function_variable)
205 return coeff*self.__class__(*([Mul(*rest)] + list(self.args[1:])))
207 @property
208 def as_integral(self):
209 return self._as_integral(self.function, self.function_variable,
210 self.transform_variable)
212 def _eval_rewrite_as_Integral(self, *args, **kwargs):
213 return self.as_integral
216def _simplify(expr, doit):
217 if doit:
218 from sympy.simplify import simplify
219 from sympy.simplify.powsimp import powdenest
220 return simplify(powdenest(piecewise_fold(expr), polar=True))
221 return expr
224def _noconds_(default):
225 """
226 This is a decorator generator for dropping convergence conditions.
228 Explanation
229 ===========
231 Suppose you define a function ``transform(*args)`` which returns a tuple of
232 the form ``(result, cond1, cond2, ...)``.
234 Decorating it ``@_noconds_(default)`` will add a new keyword argument
235 ``noconds`` to it. If ``noconds=True``, the return value will be altered to
236 be only ``result``, whereas if ``noconds=False`` the return value will not
237 be altered.
239 The default value of the ``noconds`` keyword will be ``default`` (i.e. the
240 argument of this function).
241 """
242 def make_wrapper(func):
243 @wraps(func)
244 def wrapper(*args, noconds=default, **kwargs):
245 res = func(*args, **kwargs)
246 if noconds:
247 return res[0]
248 return res
249 return wrapper
250 return make_wrapper
251_noconds = _noconds_(False)
254##########################################################################
255# Mellin Transform
256##########################################################################
258def _default_integrator(f, x):
259 return integrate(f, (x, S.Zero, S.Infinity))
262@_noconds
263def _mellin_transform(f, x, s_, integrator=_default_integrator, simplify=True):
264 """ Backend function to compute Mellin transforms. """
265 # We use a fresh dummy, because assumptions on s might drop conditions on
266 # convergence of the integral.
267 s = _dummy('s', 'mellin-transform', f)
268 F = integrator(x**(s - 1) * f, x)
270 if not F.has(Integral):
271 return _simplify(F.subs(s, s_), simplify), (S.NegativeInfinity, S.Infinity), S.true
273 if not F.is_Piecewise: # XXX can this work if integration gives continuous result now?
274 raise IntegralTransformError('Mellin', f, 'could not compute integral')
276 F, cond = F.args[0]
277 if F.has(Integral):
278 raise IntegralTransformError(
279 'Mellin', f, 'integral in unexpected form')
281 def process_conds(cond):
282 """
283 Turn ``cond`` into a strip (a, b), and auxiliary conditions.
284 """
285 from sympy.solvers.inequalities import _solve_inequality
286 a = S.NegativeInfinity
287 b = S.Infinity
288 aux = S.true
289 conds = conjuncts(to_cnf(cond))
290 t = Dummy('t', real=True)
291 for c in conds:
292 a_ = S.Infinity
293 b_ = S.NegativeInfinity
294 aux_ = []
295 for d in disjuncts(c):
296 d_ = d.replace(
297 re, lambda x: x.as_real_imag()[0]).subs(re(s), t)
298 if not d.is_Relational or \
299 d.rel_op in ('==', '!=') \
300 or d_.has(s) or not d_.has(t):
301 aux_ += [d]
302 continue
303 soln = _solve_inequality(d_, t)
304 if not soln.is_Relational or \
305 soln.rel_op in ('==', '!='):
306 aux_ += [d]
307 continue
308 if soln.lts == t:
309 b_ = Max(soln.gts, b_)
310 else:
311 a_ = Min(soln.lts, a_)
312 if a_ is not S.Infinity and a_ != b:
313 a = Max(a_, a)
314 elif b_ is not S.NegativeInfinity and b_ != a:
315 b = Min(b_, b)
316 else:
317 aux = And(aux, Or(*aux_))
318 return a, b, aux
320 conds = [process_conds(c) for c in disjuncts(cond)]
321 conds = [x for x in conds if x[2] != False]
322 conds.sort(key=lambda x: (x[0] - x[1], count_ops(x[2])))
324 if not conds:
325 raise IntegralTransformError('Mellin', f, 'no convergence found')
327 a, b, aux = conds[0]
328 return _simplify(F.subs(s, s_), simplify), (a, b), aux
331class MellinTransform(IntegralTransform):
332 """
333 Class representing unevaluated Mellin transforms.
335 For usage of this class, see the :class:`IntegralTransform` docstring.
337 For how to compute Mellin transforms, see the :func:`mellin_transform`
338 docstring.
339 """
341 _name = 'Mellin'
343 def _compute_transform(self, f, x, s, **hints):
344 return _mellin_transform(f, x, s, **hints)
346 def _as_integral(self, f, x, s):
347 return Integral(f*x**(s - 1), (x, S.Zero, S.Infinity))
349 def _collapse_extra(self, extra):
350 a = []
351 b = []
352 cond = []
353 for (sa, sb), c in extra:
354 a += [sa]
355 b += [sb]
356 cond += [c]
357 res = (Max(*a), Min(*b)), And(*cond)
358 if (res[0][0] >= res[0][1]) == True or res[1] == False:
359 raise IntegralTransformError(
360 'Mellin', None, 'no combined convergence.')
361 return res
364def mellin_transform(f, x, s, **hints):
365 r"""
366 Compute the Mellin transform `F(s)` of `f(x)`,
368 .. math :: F(s) = \int_0^\infty x^{s-1} f(x) \mathrm{d}x.
370 For all "sensible" functions, this converges absolutely in a strip
371 `a < \operatorname{Re}(s) < b`.
373 Explanation
374 ===========
376 The Mellin transform is related via change of variables to the Fourier
377 transform, and also to the (bilateral) Laplace transform.
379 This function returns ``(F, (a, b), cond)``
380 where ``F`` is the Mellin transform of ``f``, ``(a, b)`` is the fundamental strip
381 (as above), and ``cond`` are auxiliary convergence conditions.
383 If the integral cannot be computed in closed form, this function returns
384 an unevaluated :class:`MellinTransform` object.
386 For a description of possible hints, refer to the docstring of
387 :func:`sympy.integrals.transforms.IntegralTransform.doit`. If ``noconds=False``,
388 then only `F` will be returned (i.e. not ``cond``, and also not the strip
389 ``(a, b)``).
391 Examples
392 ========
394 >>> from sympy import mellin_transform, exp
395 >>> from sympy.abc import x, s
396 >>> mellin_transform(exp(-x), x, s)
397 (gamma(s), (0, oo), True)
399 See Also
400 ========
402 inverse_mellin_transform, laplace_transform, fourier_transform
403 hankel_transform, inverse_hankel_transform
404 """
405 return MellinTransform(f, x, s).doit(**hints)
408def _rewrite_sin(m_n, s, a, b):
409 """
410 Re-write the sine function ``sin(m*s + n)`` as gamma functions, compatible
411 with the strip (a, b).
413 Return ``(gamma1, gamma2, fac)`` so that ``f == fac/(gamma1 * gamma2)``.
415 Examples
416 ========
418 >>> from sympy.integrals.transforms import _rewrite_sin
419 >>> from sympy import pi, S
420 >>> from sympy.abc import s
421 >>> _rewrite_sin((pi, 0), s, 0, 1)
422 (gamma(s), gamma(1 - s), pi)
423 >>> _rewrite_sin((pi, 0), s, 1, 0)
424 (gamma(s - 1), gamma(2 - s), -pi)
425 >>> _rewrite_sin((pi, 0), s, -1, 0)
426 (gamma(s + 1), gamma(-s), -pi)
427 >>> _rewrite_sin((pi, pi/2), s, S(1)/2, S(3)/2)
428 (gamma(s - 1/2), gamma(3/2 - s), -pi)
429 >>> _rewrite_sin((pi, pi), s, 0, 1)
430 (gamma(s), gamma(1 - s), -pi)
431 >>> _rewrite_sin((2*pi, 0), s, 0, S(1)/2)
432 (gamma(2*s), gamma(1 - 2*s), pi)
433 >>> _rewrite_sin((2*pi, 0), s, S(1)/2, 1)
434 (gamma(2*s - 1), gamma(2 - 2*s), -pi)
435 """
436 # (This is a separate function because it is moderately complicated,
437 # and I want to doctest it.)
438 # We want to use pi/sin(pi*x) = gamma(x)*gamma(1-x).
439 # But there is one comlication: the gamma functions determine the
440 # inegration contour in the definition of the G-function. Usually
441 # it would not matter if this is slightly shifted, unless this way
442 # we create an undefined function!
443 # So we try to write this in such a way that the gammas are
444 # eminently on the right side of the strip.
445 m, n = m_n
447 m = expand_mul(m/pi)
448 n = expand_mul(n/pi)
449 r = ceiling(-m*a - n.as_real_imag()[0]) # Don't use re(n), does not expand
450 return gamma(m*s + n + r), gamma(1 - n - r - m*s), (-1)**r*pi
453class MellinTransformStripError(ValueError):
454 """
455 Exception raised by _rewrite_gamma. Mainly for internal use.
456 """
457 pass
460def _rewrite_gamma(f, s, a, b):
461 """
462 Try to rewrite the product f(s) as a product of gamma functions,
463 so that the inverse Mellin transform of f can be expressed as a meijer
464 G function.
466 Explanation
467 ===========
469 Return (an, ap), (bm, bq), arg, exp, fac such that
470 G((an, ap), (bm, bq), arg/z**exp)*fac is the inverse Mellin transform of f(s).
472 Raises IntegralTransformError or MellinTransformStripError on failure.
474 It is asserted that f has no poles in the fundamental strip designated by
475 (a, b). One of a and b is allowed to be None. The fundamental strip is
476 important, because it determines the inversion contour.
478 This function can handle exponentials, linear factors, trigonometric
479 functions.
481 This is a helper function for inverse_mellin_transform that will not
482 attempt any transformations on f.
484 Examples
485 ========
487 >>> from sympy.integrals.transforms import _rewrite_gamma
488 >>> from sympy.abc import s
489 >>> from sympy import oo
490 >>> _rewrite_gamma(s*(s+3)*(s-1), s, -oo, oo)
491 (([], [-3, 0, 1]), ([-2, 1, 2], []), 1, 1, -1)
492 >>> _rewrite_gamma((s-1)**2, s, -oo, oo)
493 (([], [1, 1]), ([2, 2], []), 1, 1, 1)
495 Importance of the fundamental strip:
497 >>> _rewrite_gamma(1/s, s, 0, oo)
498 (([1], []), ([], [0]), 1, 1, 1)
499 >>> _rewrite_gamma(1/s, s, None, oo)
500 (([1], []), ([], [0]), 1, 1, 1)
501 >>> _rewrite_gamma(1/s, s, 0, None)
502 (([1], []), ([], [0]), 1, 1, 1)
503 >>> _rewrite_gamma(1/s, s, -oo, 0)
504 (([], [1]), ([0], []), 1, 1, -1)
505 >>> _rewrite_gamma(1/s, s, None, 0)
506 (([], [1]), ([0], []), 1, 1, -1)
507 >>> _rewrite_gamma(1/s, s, -oo, None)
508 (([], [1]), ([0], []), 1, 1, -1)
510 >>> _rewrite_gamma(2**(-s+3), s, -oo, oo)
511 (([], []), ([], []), 1/2, 1, 8)
512 """
513 # Our strategy will be as follows:
514 # 1) Guess a constant c such that the inversion integral should be
515 # performed wrt s'=c*s (instead of plain s). Write s for s'.
516 # 2) Process all factors, rewrite them independently as gamma functions in
517 # argument s, or exponentials of s.
518 # 3) Try to transform all gamma functions s.t. they have argument
519 # a+s or a-s.
520 # 4) Check that the resulting G function parameters are valid.
521 # 5) Combine all the exponentials.
523 a_, b_ = S([a, b])
525 def left(c, is_numer):
526 """
527 Decide whether pole at c lies to the left of the fundamental strip.
528 """
529 # heuristically, this is the best chance for us to solve the inequalities
530 c = expand(re(c))
531 if a_ is None and b_ is S.Infinity:
532 return True
533 if a_ is None:
534 return c < b_
535 if b_ is None:
536 return c <= a_
537 if (c >= b_) == True:
538 return False
539 if (c <= a_) == True:
540 return True
541 if is_numer:
542 return None
543 if a_.free_symbols or b_.free_symbols or c.free_symbols:
544 return None # XXX
545 #raise IntegralTransformError('Inverse Mellin', f,
546 # 'Could not determine position of singularity %s'
547 # ' relative to fundamental strip' % c)
548 raise MellinTransformStripError('Pole inside critical strip?')
550 # 1)
551 s_multipliers = []
552 for g in f.atoms(gamma):
553 if not g.has(s):
554 continue
555 arg = g.args[0]
556 if arg.is_Add:
557 arg = arg.as_independent(s)[1]
558 coeff, _ = arg.as_coeff_mul(s)
559 s_multipliers += [coeff]
560 for g in f.atoms(sin, cos, tan, cot):
561 if not g.has(s):
562 continue
563 arg = g.args[0]
564 if arg.is_Add:
565 arg = arg.as_independent(s)[1]
566 coeff, _ = arg.as_coeff_mul(s)
567 s_multipliers += [coeff/pi]
568 s_multipliers = [Abs(x) if x.is_extended_real else x for x in s_multipliers]
569 common_coefficient = S.One
570 for x in s_multipliers:
571 if not x.is_Rational:
572 common_coefficient = x
573 break
574 s_multipliers = [x/common_coefficient for x in s_multipliers]
575 if not (all(x.is_Rational for x in s_multipliers) and
576 common_coefficient.is_extended_real):
577 raise IntegralTransformError("Gamma", None, "Nonrational multiplier")
578 s_multiplier = common_coefficient/reduce(ilcm, [S(x.q)
579 for x in s_multipliers], S.One)
580 if s_multiplier == common_coefficient:
581 if len(s_multipliers) == 0:
582 s_multiplier = common_coefficient
583 else:
584 s_multiplier = common_coefficient \
585 *reduce(igcd, [S(x.p) for x in s_multipliers])
587 f = f.subs(s, s/s_multiplier)
588 fac = S.One/s_multiplier
589 exponent = S.One/s_multiplier
590 if a_ is not None:
591 a_ *= s_multiplier
592 if b_ is not None:
593 b_ *= s_multiplier
595 # 2)
596 numer, denom = f.as_numer_denom()
597 numer = Mul.make_args(numer)
598 denom = Mul.make_args(denom)
599 args = list(zip(numer, repeat(True))) + list(zip(denom, repeat(False)))
601 facs = []
602 dfacs = []
603 # *_gammas will contain pairs (a, c) representing Gamma(a*s + c)
604 numer_gammas = []
605 denom_gammas = []
606 # exponentials will contain bases for exponentials of s
607 exponentials = []
609 def exception(fact):
610 return IntegralTransformError("Inverse Mellin", f, "Unrecognised form '%s'." % fact)
611 while args:
612 fact, is_numer = args.pop()
613 if is_numer:
614 ugammas, lgammas = numer_gammas, denom_gammas
615 ufacs = facs
616 else:
617 ugammas, lgammas = denom_gammas, numer_gammas
618 ufacs = dfacs
620 def linear_arg(arg):
621 """ Test if arg is of form a*s+b, raise exception if not. """
622 if not arg.is_polynomial(s):
623 raise exception(fact)
624 p = Poly(arg, s)
625 if p.degree() != 1:
626 raise exception(fact)
627 return p.all_coeffs()
629 # constants
630 if not fact.has(s):
631 ufacs += [fact]
632 # exponentials
633 elif fact.is_Pow or isinstance(fact, exp):
634 if fact.is_Pow:
635 base = fact.base
636 exp_ = fact.exp
637 else:
638 base = exp_polar(1)
639 exp_ = fact.exp
640 if exp_.is_Integer:
641 cond = is_numer
642 if exp_ < 0:
643 cond = not cond
644 args += [(base, cond)]*Abs(exp_)
645 continue
646 elif not base.has(s):
647 a, b = linear_arg(exp_)
648 if not is_numer:
649 base = 1/base
650 exponentials += [base**a]
651 facs += [base**b]
652 else:
653 raise exception(fact)
654 # linear factors
655 elif fact.is_polynomial(s):
656 p = Poly(fact, s)
657 if p.degree() != 1:
658 # We completely factor the poly. For this we need the roots.
659 # Now roots() only works in some cases (low degree), and CRootOf
660 # only works without parameters. So try both...
661 coeff = p.LT()[1]
662 rs = roots(p, s)
663 if len(rs) != p.degree():
664 rs = CRootOf.all_roots(p)
665 ufacs += [coeff]
666 args += [(s - c, is_numer) for c in rs]
667 continue
668 a, c = p.all_coeffs()
669 ufacs += [a]
670 c /= -a
671 # Now need to convert s - c
672 if left(c, is_numer):
673 ugammas += [(S.One, -c + 1)]
674 lgammas += [(S.One, -c)]
675 else:
676 ufacs += [-1]
677 ugammas += [(S.NegativeOne, c + 1)]
678 lgammas += [(S.NegativeOne, c)]
679 elif isinstance(fact, gamma):
680 a, b = linear_arg(fact.args[0])
681 if is_numer:
682 if (a > 0 and (left(-b/a, is_numer) == False)) or \
683 (a < 0 and (left(-b/a, is_numer) == True)):
684 raise NotImplementedError(
685 'Gammas partially over the strip.')
686 ugammas += [(a, b)]
687 elif isinstance(fact, sin):
688 # We try to re-write all trigs as gammas. This is not in
689 # general the best strategy, since sometimes this is impossible,
690 # but rewriting as exponentials would work. However trig functions
691 # in inverse mellin transforms usually all come from simplifying
692 # gamma terms, so this should work.
693 a = fact.args[0]
694 if is_numer:
695 # No problem with the poles.
696 gamma1, gamma2, fac_ = gamma(a/pi), gamma(1 - a/pi), pi
697 else:
698 gamma1, gamma2, fac_ = _rewrite_sin(linear_arg(a), s, a_, b_)
699 args += [(gamma1, not is_numer), (gamma2, not is_numer)]
700 ufacs += [fac_]
701 elif isinstance(fact, tan):
702 a = fact.args[0]
703 args += [(sin(a, evaluate=False), is_numer),
704 (sin(pi/2 - a, evaluate=False), not is_numer)]
705 elif isinstance(fact, cos):
706 a = fact.args[0]
707 args += [(sin(pi/2 - a, evaluate=False), is_numer)]
708 elif isinstance(fact, cot):
709 a = fact.args[0]
710 args += [(sin(pi/2 - a, evaluate=False), is_numer),
711 (sin(a, evaluate=False), not is_numer)]
712 else:
713 raise exception(fact)
715 fac *= Mul(*facs)/Mul(*dfacs)
717 # 3)
718 an, ap, bm, bq = [], [], [], []
719 for gammas, plus, minus, is_numer in [(numer_gammas, an, bm, True),
720 (denom_gammas, bq, ap, False)]:
721 while gammas:
722 a, c = gammas.pop()
723 if a != -1 and a != +1:
724 # We use the gamma function multiplication theorem.
725 p = Abs(S(a))
726 newa = a/p
727 newc = c/p
728 if not a.is_Integer:
729 raise TypeError("a is not an integer")
730 for k in range(p):
731 gammas += [(newa, newc + k/p)]
732 if is_numer:
733 fac *= (2*pi)**((1 - p)/2) * p**(c - S.Half)
734 exponentials += [p**a]
735 else:
736 fac /= (2*pi)**((1 - p)/2) * p**(c - S.Half)
737 exponentials += [p**(-a)]
738 continue
739 if a == +1:
740 plus.append(1 - c)
741 else:
742 minus.append(c)
744 # 4)
745 # TODO
747 # 5)
748 arg = Mul(*exponentials)
750 # for testability, sort the arguments
751 an.sort(key=default_sort_key)
752 ap.sort(key=default_sort_key)
753 bm.sort(key=default_sort_key)
754 bq.sort(key=default_sort_key)
756 return (an, ap), (bm, bq), arg, exponent, fac
759@_noconds_(True)
760def _inverse_mellin_transform(F, s, x_, strip, as_meijerg=False):
761 """ A helper for the real inverse_mellin_transform function, this one here
762 assumes x to be real and positive. """
763 x = _dummy('t', 'inverse-mellin-transform', F, positive=True)
764 # Actually, we won't try integration at all. Instead we use the definition
765 # of the Meijer G function as a fairly general inverse mellin transform.
766 F = F.rewrite(gamma)
767 for g in [factor(F), expand_mul(F), expand(F)]:
768 if g.is_Add:
769 # do all terms separately
770 ress = [_inverse_mellin_transform(G, s, x, strip, as_meijerg,
771 noconds=False)
772 for G in g.args]
773 conds = [p[1] for p in ress]
774 ress = [p[0] for p in ress]
775 res = Add(*ress)
776 if not as_meijerg:
777 res = factor(res, gens=res.atoms(Heaviside))
778 return res.subs(x, x_), And(*conds)
780 try:
781 a, b, C, e, fac = _rewrite_gamma(g, s, strip[0], strip[1])
782 except IntegralTransformError:
783 continue
784 try:
785 G = meijerg(a, b, C/x**e)
786 except ValueError:
787 continue
788 if as_meijerg:
789 h = G
790 else:
791 try:
792 from sympy.simplify import hyperexpand
793 h = hyperexpand(G)
794 except NotImplementedError:
795 raise IntegralTransformError(
796 'Inverse Mellin', F, 'Could not calculate integral')
798 if h.is_Piecewise and len(h.args) == 3:
799 # XXX we break modularity here!
800 h = Heaviside(x - Abs(C))*h.args[0].args[0] \
801 + Heaviside(Abs(C) - x)*h.args[1].args[0]
802 # We must ensure that the integral along the line we want converges,
803 # and return that value.
804 # See [L], 5.2
805 cond = [Abs(arg(G.argument)) < G.delta*pi]
806 # Note: we allow ">=" here, this corresponds to convergence if we let
807 # limits go to oo symmetrically. ">" corresponds to absolute convergence.
808 cond += [And(Or(len(G.ap) != len(G.bq), 0 >= re(G.nu) + 1),
809 Abs(arg(G.argument)) == G.delta*pi)]
810 cond = Or(*cond)
811 if cond == False:
812 raise IntegralTransformError(
813 'Inverse Mellin', F, 'does not converge')
814 return (h*fac).subs(x, x_), cond
816 raise IntegralTransformError('Inverse Mellin', F, '')
818_allowed = None
821class InverseMellinTransform(IntegralTransform):
822 """
823 Class representing unevaluated inverse Mellin transforms.
825 For usage of this class, see the :class:`IntegralTransform` docstring.
827 For how to compute inverse Mellin transforms, see the
828 :func:`inverse_mellin_transform` docstring.
829 """
831 _name = 'Inverse Mellin'
832 _none_sentinel = Dummy('None')
833 _c = Dummy('c')
835 def __new__(cls, F, s, x, a, b, **opts):
836 if a is None:
837 a = InverseMellinTransform._none_sentinel
838 if b is None:
839 b = InverseMellinTransform._none_sentinel
840 return IntegralTransform.__new__(cls, F, s, x, a, b, **opts)
842 @property
843 def fundamental_strip(self):
844 a, b = self.args[3], self.args[4]
845 if a is InverseMellinTransform._none_sentinel:
846 a = None
847 if b is InverseMellinTransform._none_sentinel:
848 b = None
849 return a, b
851 def _compute_transform(self, F, s, x, **hints):
852 # IntegralTransform's doit will cause this hint to exist, but
853 # InverseMellinTransform should ignore it
854 hints.pop('simplify', True)
855 global _allowed
856 if _allowed is None:
857 _allowed = {
858 exp, gamma, sin, cos, tan, cot, cosh, sinh, tanh, coth,
859 factorial, rf}
860 for f in postorder_traversal(F):
861 if f.is_Function and f.has(s) and f.func not in _allowed:
862 raise IntegralTransformError('Inverse Mellin', F,
863 'Component %s not recognised.' % f)
864 strip = self.fundamental_strip
865 return _inverse_mellin_transform(F, s, x, strip, **hints)
867 def _as_integral(self, F, s, x):
868 c = self.__class__._c
869 return Integral(F*x**(-s), (s, c - S.ImaginaryUnit*S.Infinity, c +
870 S.ImaginaryUnit*S.Infinity))/(2*S.Pi*S.ImaginaryUnit)
873def inverse_mellin_transform(F, s, x, strip, **hints):
874 r"""
875 Compute the inverse Mellin transform of `F(s)` over the fundamental
876 strip given by ``strip=(a, b)``.
878 Explanation
879 ===========
881 This can be defined as
883 .. math:: f(x) = \frac{1}{2\pi i} \int_{c - i\infty}^{c + i\infty} x^{-s} F(s) \mathrm{d}s,
885 for any `c` in the fundamental strip. Under certain regularity
886 conditions on `F` and/or `f`,
887 this recovers `f` from its Mellin transform `F`
888 (and vice versa), for positive real `x`.
890 One of `a` or `b` may be passed as ``None``; a suitable `c` will be
891 inferred.
893 If the integral cannot be computed in closed form, this function returns
894 an unevaluated :class:`InverseMellinTransform` object.
896 Note that this function will assume x to be positive and real, regardless
897 of the SymPy assumptions!
899 For a description of possible hints, refer to the docstring of
900 :func:`sympy.integrals.transforms.IntegralTransform.doit`.
902 Examples
903 ========
905 >>> from sympy import inverse_mellin_transform, oo, gamma
906 >>> from sympy.abc import x, s
907 >>> inverse_mellin_transform(gamma(s), s, x, (0, oo))
908 exp(-x)
910 The fundamental strip matters:
912 >>> f = 1/(s**2 - 1)
913 >>> inverse_mellin_transform(f, s, x, (-oo, -1))
914 x*(1 - 1/x**2)*Heaviside(x - 1)/2
915 >>> inverse_mellin_transform(f, s, x, (-1, 1))
916 -x*Heaviside(1 - x)/2 - Heaviside(x - 1)/(2*x)
917 >>> inverse_mellin_transform(f, s, x, (1, oo))
918 (1/2 - x**2/2)*Heaviside(1 - x)/x
920 See Also
921 ========
923 mellin_transform
924 hankel_transform, inverse_hankel_transform
925 """
926 return InverseMellinTransform(F, s, x, strip[0], strip[1]).doit(**hints)
929##########################################################################
930# Fourier Transform
931##########################################################################
933@_noconds_(True)
934def _fourier_transform(f, x, k, a, b, name, simplify=True):
935 r"""
936 Compute a general Fourier-type transform
938 .. math::
940 F(k) = a \int_{-\infty}^{\infty} e^{bixk} f(x)\, dx.
942 For suitable choice of *a* and *b*, this reduces to the standard Fourier
943 and inverse Fourier transforms.
944 """
945 F = integrate(a*f*exp(b*S.ImaginaryUnit*x*k), (x, S.NegativeInfinity, S.Infinity))
947 if not F.has(Integral):
948 return _simplify(F, simplify), S.true
950 integral_f = integrate(f, (x, S.NegativeInfinity, S.Infinity))
951 if integral_f in (S.NegativeInfinity, S.Infinity, S.NaN) or integral_f.has(Integral):
952 raise IntegralTransformError(name, f, 'function not integrable on real axis')
954 if not F.is_Piecewise:
955 raise IntegralTransformError(name, f, 'could not compute integral')
957 F, cond = F.args[0]
958 if F.has(Integral):
959 raise IntegralTransformError(name, f, 'integral in unexpected form')
961 return _simplify(F, simplify), cond
964class FourierTypeTransform(IntegralTransform):
965 """ Base class for Fourier transforms."""
967 def a(self):
968 raise NotImplementedError(
969 "Class %s must implement a(self) but does not" % self.__class__)
971 def b(self):
972 raise NotImplementedError(
973 "Class %s must implement b(self) but does not" % self.__class__)
975 def _compute_transform(self, f, x, k, **hints):
976 return _fourier_transform(f, x, k,
977 self.a(), self.b(),
978 self.__class__._name, **hints)
980 def _as_integral(self, f, x, k):
981 a = self.a()
982 b = self.b()
983 return Integral(a*f*exp(b*S.ImaginaryUnit*x*k), (x, S.NegativeInfinity, S.Infinity))
986class FourierTransform(FourierTypeTransform):
987 """
988 Class representing unevaluated Fourier transforms.
990 For usage of this class, see the :class:`IntegralTransform` docstring.
992 For how to compute Fourier transforms, see the :func:`fourier_transform`
993 docstring.
994 """
996 _name = 'Fourier'
998 def a(self):
999 return 1
1001 def b(self):
1002 return -2*S.Pi
1005def fourier_transform(f, x, k, **hints):
1006 r"""
1007 Compute the unitary, ordinary-frequency Fourier transform of ``f``, defined
1008 as
1010 .. math:: F(k) = \int_{-\infty}^\infty f(x) e^{-2\pi i x k} \mathrm{d} x.
1012 Explanation
1013 ===========
1015 If the transform cannot be computed in closed form, this
1016 function returns an unevaluated :class:`FourierTransform` object.
1018 For other Fourier transform conventions, see the function
1019 :func:`sympy.integrals.transforms._fourier_transform`.
1021 For a description of possible hints, refer to the docstring of
1022 :func:`sympy.integrals.transforms.IntegralTransform.doit`.
1023 Note that for this transform, by default ``noconds=True``.
1025 Examples
1026 ========
1028 >>> from sympy import fourier_transform, exp
1029 >>> from sympy.abc import x, k
1030 >>> fourier_transform(exp(-x**2), x, k)
1031 sqrt(pi)*exp(-pi**2*k**2)
1032 >>> fourier_transform(exp(-x**2), x, k, noconds=False)
1033 (sqrt(pi)*exp(-pi**2*k**2), True)
1035 See Also
1036 ========
1038 inverse_fourier_transform
1039 sine_transform, inverse_sine_transform
1040 cosine_transform, inverse_cosine_transform
1041 hankel_transform, inverse_hankel_transform
1042 mellin_transform, laplace_transform
1043 """
1044 return FourierTransform(f, x, k).doit(**hints)
1047class InverseFourierTransform(FourierTypeTransform):
1048 """
1049 Class representing unevaluated inverse Fourier transforms.
1051 For usage of this class, see the :class:`IntegralTransform` docstring.
1053 For how to compute inverse Fourier transforms, see the
1054 :func:`inverse_fourier_transform` docstring.
1055 """
1057 _name = 'Inverse Fourier'
1059 def a(self):
1060 return 1
1062 def b(self):
1063 return 2*S.Pi
1066def inverse_fourier_transform(F, k, x, **hints):
1067 r"""
1068 Compute the unitary, ordinary-frequency inverse Fourier transform of `F`,
1069 defined as
1071 .. math:: f(x) = \int_{-\infty}^\infty F(k) e^{2\pi i x k} \mathrm{d} k.
1073 Explanation
1074 ===========
1076 If the transform cannot be computed in closed form, this
1077 function returns an unevaluated :class:`InverseFourierTransform` object.
1079 For other Fourier transform conventions, see the function
1080 :func:`sympy.integrals.transforms._fourier_transform`.
1082 For a description of possible hints, refer to the docstring of
1083 :func:`sympy.integrals.transforms.IntegralTransform.doit`.
1084 Note that for this transform, by default ``noconds=True``.
1086 Examples
1087 ========
1089 >>> from sympy import inverse_fourier_transform, exp, sqrt, pi
1090 >>> from sympy.abc import x, k
1091 >>> inverse_fourier_transform(sqrt(pi)*exp(-(pi*k)**2), k, x)
1092 exp(-x**2)
1093 >>> inverse_fourier_transform(sqrt(pi)*exp(-(pi*k)**2), k, x, noconds=False)
1094 (exp(-x**2), True)
1096 See Also
1097 ========
1099 fourier_transform
1100 sine_transform, inverse_sine_transform
1101 cosine_transform, inverse_cosine_transform
1102 hankel_transform, inverse_hankel_transform
1103 mellin_transform, laplace_transform
1104 """
1105 return InverseFourierTransform(F, k, x).doit(**hints)
1108##########################################################################
1109# Fourier Sine and Cosine Transform
1110##########################################################################
1112@_noconds_(True)
1113def _sine_cosine_transform(f, x, k, a, b, K, name, simplify=True):
1114 """
1115 Compute a general sine or cosine-type transform
1116 F(k) = a int_0^oo b*sin(x*k) f(x) dx.
1117 F(k) = a int_0^oo b*cos(x*k) f(x) dx.
1119 For suitable choice of a and b, this reduces to the standard sine/cosine
1120 and inverse sine/cosine transforms.
1121 """
1122 F = integrate(a*f*K(b*x*k), (x, S.Zero, S.Infinity))
1124 if not F.has(Integral):
1125 return _simplify(F, simplify), S.true
1127 if not F.is_Piecewise:
1128 raise IntegralTransformError(name, f, 'could not compute integral')
1130 F, cond = F.args[0]
1131 if F.has(Integral):
1132 raise IntegralTransformError(name, f, 'integral in unexpected form')
1134 return _simplify(F, simplify), cond
1137class SineCosineTypeTransform(IntegralTransform):
1138 """
1139 Base class for sine and cosine transforms.
1140 Specify cls._kern.
1141 """
1143 def a(self):
1144 raise NotImplementedError(
1145 "Class %s must implement a(self) but does not" % self.__class__)
1147 def b(self):
1148 raise NotImplementedError(
1149 "Class %s must implement b(self) but does not" % self.__class__)
1152 def _compute_transform(self, f, x, k, **hints):
1153 return _sine_cosine_transform(f, x, k,
1154 self.a(), self.b(),
1155 self.__class__._kern,
1156 self.__class__._name, **hints)
1158 def _as_integral(self, f, x, k):
1159 a = self.a()
1160 b = self.b()
1161 K = self.__class__._kern
1162 return Integral(a*f*K(b*x*k), (x, S.Zero, S.Infinity))
1165class SineTransform(SineCosineTypeTransform):
1166 """
1167 Class representing unevaluated sine transforms.
1169 For usage of this class, see the :class:`IntegralTransform` docstring.
1171 For how to compute sine transforms, see the :func:`sine_transform`
1172 docstring.
1173 """
1175 _name = 'Sine'
1176 _kern = sin
1178 def a(self):
1179 return sqrt(2)/sqrt(pi)
1181 def b(self):
1182 return S.One
1185def sine_transform(f, x, k, **hints):
1186 r"""
1187 Compute the unitary, ordinary-frequency sine transform of `f`, defined
1188 as
1190 .. math:: F(k) = \sqrt{\frac{2}{\pi}} \int_{0}^\infty f(x) \sin(2\pi x k) \mathrm{d} x.
1192 Explanation
1193 ===========
1195 If the transform cannot be computed in closed form, this
1196 function returns an unevaluated :class:`SineTransform` object.
1198 For a description of possible hints, refer to the docstring of
1199 :func:`sympy.integrals.transforms.IntegralTransform.doit`.
1200 Note that for this transform, by default ``noconds=True``.
1202 Examples
1203 ========
1205 >>> from sympy import sine_transform, exp
1206 >>> from sympy.abc import x, k, a
1207 >>> sine_transform(x*exp(-a*x**2), x, k)
1208 sqrt(2)*k*exp(-k**2/(4*a))/(4*a**(3/2))
1209 >>> sine_transform(x**(-a), x, k)
1210 2**(1/2 - a)*k**(a - 1)*gamma(1 - a/2)/gamma(a/2 + 1/2)
1212 See Also
1213 ========
1215 fourier_transform, inverse_fourier_transform
1216 inverse_sine_transform
1217 cosine_transform, inverse_cosine_transform
1218 hankel_transform, inverse_hankel_transform
1219 mellin_transform, laplace_transform
1220 """
1221 return SineTransform(f, x, k).doit(**hints)
1224class InverseSineTransform(SineCosineTypeTransform):
1225 """
1226 Class representing unevaluated inverse sine transforms.
1228 For usage of this class, see the :class:`IntegralTransform` docstring.
1230 For how to compute inverse sine transforms, see the
1231 :func:`inverse_sine_transform` docstring.
1232 """
1234 _name = 'Inverse Sine'
1235 _kern = sin
1237 def a(self):
1238 return sqrt(2)/sqrt(pi)
1240 def b(self):
1241 return S.One
1244def inverse_sine_transform(F, k, x, **hints):
1245 r"""
1246 Compute the unitary, ordinary-frequency inverse sine transform of `F`,
1247 defined as
1249 .. math:: f(x) = \sqrt{\frac{2}{\pi}} \int_{0}^\infty F(k) \sin(2\pi x k) \mathrm{d} k.
1251 Explanation
1252 ===========
1254 If the transform cannot be computed in closed form, this
1255 function returns an unevaluated :class:`InverseSineTransform` object.
1257 For a description of possible hints, refer to the docstring of
1258 :func:`sympy.integrals.transforms.IntegralTransform.doit`.
1259 Note that for this transform, by default ``noconds=True``.
1261 Examples
1262 ========
1264 >>> from sympy import inverse_sine_transform, exp, sqrt, gamma
1265 >>> from sympy.abc import x, k, a
1266 >>> inverse_sine_transform(2**((1-2*a)/2)*k**(a - 1)*
1267 ... gamma(-a/2 + 1)/gamma((a+1)/2), k, x)
1268 x**(-a)
1269 >>> inverse_sine_transform(sqrt(2)*k*exp(-k**2/(4*a))/(4*sqrt(a)**3), k, x)
1270 x*exp(-a*x**2)
1272 See Also
1273 ========
1275 fourier_transform, inverse_fourier_transform
1276 sine_transform
1277 cosine_transform, inverse_cosine_transform
1278 hankel_transform, inverse_hankel_transform
1279 mellin_transform, laplace_transform
1280 """
1281 return InverseSineTransform(F, k, x).doit(**hints)
1284class CosineTransform(SineCosineTypeTransform):
1285 """
1286 Class representing unevaluated cosine transforms.
1288 For usage of this class, see the :class:`IntegralTransform` docstring.
1290 For how to compute cosine transforms, see the :func:`cosine_transform`
1291 docstring.
1292 """
1294 _name = 'Cosine'
1295 _kern = cos
1297 def a(self):
1298 return sqrt(2)/sqrt(pi)
1300 def b(self):
1301 return S.One
1304def cosine_transform(f, x, k, **hints):
1305 r"""
1306 Compute the unitary, ordinary-frequency cosine transform of `f`, defined
1307 as
1309 .. math:: F(k) = \sqrt{\frac{2}{\pi}} \int_{0}^\infty f(x) \cos(2\pi x k) \mathrm{d} x.
1311 Explanation
1312 ===========
1314 If the transform cannot be computed in closed form, this
1315 function returns an unevaluated :class:`CosineTransform` object.
1317 For a description of possible hints, refer to the docstring of
1318 :func:`sympy.integrals.transforms.IntegralTransform.doit`.
1319 Note that for this transform, by default ``noconds=True``.
1321 Examples
1322 ========
1324 >>> from sympy import cosine_transform, exp, sqrt, cos
1325 >>> from sympy.abc import x, k, a
1326 >>> cosine_transform(exp(-a*x), x, k)
1327 sqrt(2)*a/(sqrt(pi)*(a**2 + k**2))
1328 >>> cosine_transform(exp(-a*sqrt(x))*cos(a*sqrt(x)), x, k)
1329 a*exp(-a**2/(2*k))/(2*k**(3/2))
1331 See Also
1332 ========
1334 fourier_transform, inverse_fourier_transform,
1335 sine_transform, inverse_sine_transform
1336 inverse_cosine_transform
1337 hankel_transform, inverse_hankel_transform
1338 mellin_transform, laplace_transform
1339 """
1340 return CosineTransform(f, x, k).doit(**hints)
1343class InverseCosineTransform(SineCosineTypeTransform):
1344 """
1345 Class representing unevaluated inverse cosine transforms.
1347 For usage of this class, see the :class:`IntegralTransform` docstring.
1349 For how to compute inverse cosine transforms, see the
1350 :func:`inverse_cosine_transform` docstring.
1351 """
1353 _name = 'Inverse Cosine'
1354 _kern = cos
1356 def a(self):
1357 return sqrt(2)/sqrt(pi)
1359 def b(self):
1360 return S.One
1363def inverse_cosine_transform(F, k, x, **hints):
1364 r"""
1365 Compute the unitary, ordinary-frequency inverse cosine transform of `F`,
1366 defined as
1368 .. math:: f(x) = \sqrt{\frac{2}{\pi}} \int_{0}^\infty F(k) \cos(2\pi x k) \mathrm{d} k.
1370 Explanation
1371 ===========
1373 If the transform cannot be computed in closed form, this
1374 function returns an unevaluated :class:`InverseCosineTransform` object.
1376 For a description of possible hints, refer to the docstring of
1377 :func:`sympy.integrals.transforms.IntegralTransform.doit`.
1378 Note that for this transform, by default ``noconds=True``.
1380 Examples
1381 ========
1383 >>> from sympy import inverse_cosine_transform, sqrt, pi
1384 >>> from sympy.abc import x, k, a
1385 >>> inverse_cosine_transform(sqrt(2)*a/(sqrt(pi)*(a**2 + k**2)), k, x)
1386 exp(-a*x)
1387 >>> inverse_cosine_transform(1/sqrt(k), k, x)
1388 1/sqrt(x)
1390 See Also
1391 ========
1393 fourier_transform, inverse_fourier_transform,
1394 sine_transform, inverse_sine_transform
1395 cosine_transform
1396 hankel_transform, inverse_hankel_transform
1397 mellin_transform, laplace_transform
1398 """
1399 return InverseCosineTransform(F, k, x).doit(**hints)
1402##########################################################################
1403# Hankel Transform
1404##########################################################################
1406@_noconds_(True)
1407def _hankel_transform(f, r, k, nu, name, simplify=True):
1408 r"""
1409 Compute a general Hankel transform
1411 .. math:: F_\nu(k) = \int_{0}^\infty f(r) J_\nu(k r) r \mathrm{d} r.
1412 """
1413 F = integrate(f*besselj(nu, k*r)*r, (r, S.Zero, S.Infinity))
1415 if not F.has(Integral):
1416 return _simplify(F, simplify), S.true
1418 if not F.is_Piecewise:
1419 raise IntegralTransformError(name, f, 'could not compute integral')
1421 F, cond = F.args[0]
1422 if F.has(Integral):
1423 raise IntegralTransformError(name, f, 'integral in unexpected form')
1425 return _simplify(F, simplify), cond
1428class HankelTypeTransform(IntegralTransform):
1429 """
1430 Base class for Hankel transforms.
1431 """
1433 def doit(self, **hints):
1434 return self._compute_transform(self.function,
1435 self.function_variable,
1436 self.transform_variable,
1437 self.args[3],
1438 **hints)
1440 def _compute_transform(self, f, r, k, nu, **hints):
1441 return _hankel_transform(f, r, k, nu, self._name, **hints)
1443 def _as_integral(self, f, r, k, nu):
1444 return Integral(f*besselj(nu, k*r)*r, (r, S.Zero, S.Infinity))
1446 @property
1447 def as_integral(self):
1448 return self._as_integral(self.function,
1449 self.function_variable,
1450 self.transform_variable,
1451 self.args[3])
1454class HankelTransform(HankelTypeTransform):
1455 """
1456 Class representing unevaluated Hankel transforms.
1458 For usage of this class, see the :class:`IntegralTransform` docstring.
1460 For how to compute Hankel transforms, see the :func:`hankel_transform`
1461 docstring.
1462 """
1464 _name = 'Hankel'
1467def hankel_transform(f, r, k, nu, **hints):
1468 r"""
1469 Compute the Hankel transform of `f`, defined as
1471 .. math:: F_\nu(k) = \int_{0}^\infty f(r) J_\nu(k r) r \mathrm{d} r.
1473 Explanation
1474 ===========
1476 If the transform cannot be computed in closed form, this
1477 function returns an unevaluated :class:`HankelTransform` object.
1479 For a description of possible hints, refer to the docstring of
1480 :func:`sympy.integrals.transforms.IntegralTransform.doit`.
1481 Note that for this transform, by default ``noconds=True``.
1483 Examples
1484 ========
1486 >>> from sympy import hankel_transform, inverse_hankel_transform
1487 >>> from sympy import exp
1488 >>> from sympy.abc import r, k, m, nu, a
1490 >>> ht = hankel_transform(1/r**m, r, k, nu)
1491 >>> ht
1492 2*k**(m - 2)*gamma(-m/2 + nu/2 + 1)/(2**m*gamma(m/2 + nu/2))
1494 >>> inverse_hankel_transform(ht, k, r, nu)
1495 r**(-m)
1497 >>> ht = hankel_transform(exp(-a*r), r, k, 0)
1498 >>> ht
1499 a/(k**3*(a**2/k**2 + 1)**(3/2))
1501 >>> inverse_hankel_transform(ht, k, r, 0)
1502 exp(-a*r)
1504 See Also
1505 ========
1507 fourier_transform, inverse_fourier_transform
1508 sine_transform, inverse_sine_transform
1509 cosine_transform, inverse_cosine_transform
1510 inverse_hankel_transform
1511 mellin_transform, laplace_transform
1512 """
1513 return HankelTransform(f, r, k, nu).doit(**hints)
1516class InverseHankelTransform(HankelTypeTransform):
1517 """
1518 Class representing unevaluated inverse Hankel transforms.
1520 For usage of this class, see the :class:`IntegralTransform` docstring.
1522 For how to compute inverse Hankel transforms, see the
1523 :func:`inverse_hankel_transform` docstring.
1524 """
1526 _name = 'Inverse Hankel'
1529def inverse_hankel_transform(F, k, r, nu, **hints):
1530 r"""
1531 Compute the inverse Hankel transform of `F` defined as
1533 .. math:: f(r) = \int_{0}^\infty F_\nu(k) J_\nu(k r) k \mathrm{d} k.
1535 Explanation
1536 ===========
1538 If the transform cannot be computed in closed form, this
1539 function returns an unevaluated :class:`InverseHankelTransform` object.
1541 For a description of possible hints, refer to the docstring of
1542 :func:`sympy.integrals.transforms.IntegralTransform.doit`.
1543 Note that for this transform, by default ``noconds=True``.
1545 Examples
1546 ========
1548 >>> from sympy import hankel_transform, inverse_hankel_transform
1549 >>> from sympy import exp
1550 >>> from sympy.abc import r, k, m, nu, a
1552 >>> ht = hankel_transform(1/r**m, r, k, nu)
1553 >>> ht
1554 2*k**(m - 2)*gamma(-m/2 + nu/2 + 1)/(2**m*gamma(m/2 + nu/2))
1556 >>> inverse_hankel_transform(ht, k, r, nu)
1557 r**(-m)
1559 >>> ht = hankel_transform(exp(-a*r), r, k, 0)
1560 >>> ht
1561 a/(k**3*(a**2/k**2 + 1)**(3/2))
1563 >>> inverse_hankel_transform(ht, k, r, 0)
1564 exp(-a*r)
1566 See Also
1567 ========
1569 fourier_transform, inverse_fourier_transform
1570 sine_transform, inverse_sine_transform
1571 cosine_transform, inverse_cosine_transform
1572 hankel_transform
1573 mellin_transform, laplace_transform
1574 """
1575 return InverseHankelTransform(F, k, r, nu).doit(**hints)
1578##########################################################################
1579# Laplace Transform
1580##########################################################################
1582# Stub classes and functions that used to be here
1583import sympy.integrals.laplace as _laplace
1585LaplaceTransform = _laplace.LaplaceTransform
1586laplace_transform = _laplace.laplace_transform
1587InverseLaplaceTransform = _laplace.InverseLaplaceTransform
1588inverse_laplace_transform = _laplace.inverse_laplace_transform