""" Integral Transforms """ from functools import reduce, wraps from itertools import repeat from sympy.core import S, pi, I from sympy.core.add import Add from sympy.core.function import (AppliedUndef, count_ops, Derivative, expand, expand_complex, expand_mul, Function, Lambda, WildFunction) from sympy.core.mul import Mul from sympy.core.numbers import igcd, ilcm from sympy.core.relational import _canonical, Ge, Gt, Lt, Unequality, Eq from sympy.core.sorting import default_sort_key, ordered from sympy.core.symbol import Dummy, symbols, Wild from sympy.core.traversal import postorder_traversal from sympy.functions.combinatorial.factorials import factorial, rf from sympy.functions.elementary.complexes import (re, arg, Abs, polar_lift, periodic_argument) from sympy.functions.elementary.exponential import exp, log, exp_polar from sympy.functions.elementary.hyperbolic import cosh, coth, sinh, tanh, asinh from sympy.functions.elementary.integers import ceiling from sympy.functions.elementary.miscellaneous import Max, Min, sqrt from sympy.functions.elementary.piecewise import Piecewise, piecewise_fold from sympy.functions.elementary.trigonometric import cos, cot, sin, tan, atan from sympy.functions.special.bessel import besseli, besselj, besselk, bessely from sympy.functions.special.delta_functions import DiracDelta, Heaviside from sympy.functions.special.error_functions import erf, erfc, Ei from sympy.functions.special.gamma_functions import digamma, gamma, lowergamma from sympy.functions.special.hyper import meijerg from sympy.integrals import integrate, Integral from sympy.integrals.meijerint import _dummy from sympy.logic.boolalg import to_cnf, conjuncts, disjuncts, Or, And from sympy.matrices.matrices import MatrixBase from sympy.polys.matrices.linsolve import _lin_eq2dict, PolyNonlinearError from sympy.polys.polyroots import roots from sympy.polys.polytools import factor, Poly from sympy.polys.rationaltools import together from sympy.polys.rootoftools import CRootOf, RootSum from sympy.utilities.exceptions import (sympy_deprecation_warning, SymPyDeprecationWarning, ignore_warnings) from sympy.utilities.iterables import iterable from sympy.utilities.misc import debug ########################################################################## # Helpers / Utilities ########################################################################## class IntegralTransformError(NotImplementedError): """ Exception raised in relation to problems computing transforms. Explanation =========== This class is mostly used internally; if integrals cannot be computed objects representing unevaluated transforms are usually returned. The hint ``needeval=True`` can be used to disable returning transform objects, and instead raise this exception if an integral cannot be computed. """ def __init__(self, transform, function, msg): super().__init__( "%s Transform could not be computed: %s." % (transform, msg)) self.function = function class IntegralTransform(Function): """ Base class for integral transforms. Explanation =========== This class represents unevaluated transforms. To implement a concrete transform, derive from this class and implement the ``_compute_transform(f, x, s, **hints)`` and ``_as_integral(f, x, s)`` functions. If the transform cannot be computed, raise :obj:`IntegralTransformError`. Also set ``cls._name``. For instance, >>> from sympy import LaplaceTransform >>> LaplaceTransform._name 'Laplace' Implement ``self._collapse_extra`` if your function returns more than just a number and possibly a convergence condition. """ @property def function(self): """ The function to be transformed. """ return self.args[0] @property def function_variable(self): """ The dependent variable of the function to be transformed. """ return self.args[1] @property def transform_variable(self): """ The independent transform variable. """ return self.args[2] @property def free_symbols(self): """ This method returns the symbols that will exist when the transform is evaluated. """ return self.function.free_symbols.union({self.transform_variable}) \ - {self.function_variable} def _compute_transform(self, f, x, s, **hints): raise NotImplementedError def _as_integral(self, f, x, s): raise NotImplementedError def _collapse_extra(self, extra): cond = And(*extra) if cond == False: raise IntegralTransformError(self.__class__.name, None, '') return cond def _try_directly(self, **hints): T = None try_directly = not any(func.has(self.function_variable) for func in self.function.atoms(AppliedUndef)) if try_directly: try: T = self._compute_transform(self.function, self.function_variable, self.transform_variable, **hints) except IntegralTransformError: T = None fn = self.function if not fn.is_Add: fn = expand_mul(fn) return fn, T def doit(self, **hints): """ Try to evaluate the transform in closed form. Explanation =========== This general function handles linearity, but apart from that leaves pretty much everything to _compute_transform. Standard hints are the following: - ``simplify``: whether or not to simplify the result - ``noconds``: if True, do not return convergence conditions - ``needeval``: if True, raise IntegralTransformError instead of returning IntegralTransform objects The default values of these hints depend on the concrete transform, usually the default is ``(simplify, noconds, needeval) = (True, False, False)``. """ needeval = hints.pop('needeval', False) simplify = hints.pop('simplify', True) hints['simplify'] = simplify fn, T = self._try_directly(**hints) if T is not None: return T if fn.is_Add: hints['needeval'] = needeval res = [self.__class__(*([x] + list(self.args[1:]))).doit(**hints) for x in fn.args] extra = [] ress = [] for x in res: if not isinstance(x, tuple): x = [x] ress.append(x[0]) if len(x) == 2: # only a condition extra.append(x[1]) elif len(x) > 2: # some region parameters and a condition (Mellin, Laplace) extra += [x[1:]] if simplify==True: res = Add(*ress).simplify() else: res = Add(*ress) if not extra: return res try: extra = self._collapse_extra(extra) if iterable(extra): return tuple([res]) + tuple(extra) else: return (res, extra) except IntegralTransformError: pass if needeval: raise IntegralTransformError( self.__class__._name, self.function, 'needeval') # TODO handle derivatives etc # pull out constant coefficients coeff, rest = fn.as_coeff_mul(self.function_variable) return coeff*self.__class__(*([Mul(*rest)] + list(self.args[1:]))) @property def as_integral(self): return self._as_integral(self.function, self.function_variable, self.transform_variable) def _eval_rewrite_as_Integral(self, *args, **kwargs): return self.as_integral def _simplify(expr, doit): if doit: from sympy.simplify import simplify from sympy.simplify.powsimp import powdenest return simplify(powdenest(piecewise_fold(expr), polar=True)) return expr def _noconds_(default): """ This is a decorator generator for dropping convergence conditions. Explanation =========== Suppose you define a function ``transform(*args)`` which returns a tuple of the form ``(result, cond1, cond2, ...)``. Decorating it ``@_noconds_(default)`` will add a new keyword argument ``noconds`` to it. If ``noconds=True``, the return value will be altered to be only ``result``, whereas if ``noconds=False`` the return value will not be altered. The default value of the ``noconds`` keyword will be ``default`` (i.e. the argument of this function). """ def make_wrapper(func): @wraps(func) def wrapper(*args, noconds=default, **kwargs): res = func(*args, **kwargs) if noconds: return res[0] return res return wrapper return make_wrapper _noconds = _noconds_(False) ########################################################################## # Mellin Transform ########################################################################## def _default_integrator(f, x): return integrate(f, (x, S.Zero, S.Infinity)) @_noconds def _mellin_transform(f, x, s_, integrator=_default_integrator, simplify=True): """ Backend function to compute Mellin transforms. """ # We use a fresh dummy, because assumptions on s might drop conditions on # convergence of the integral. s = _dummy('s', 'mellin-transform', f) F = integrator(x**(s - 1) * f, x) if not F.has(Integral): return _simplify(F.subs(s, s_), simplify), (S.NegativeInfinity, S.Infinity), S.true if not F.is_Piecewise: # XXX can this work if integration gives continuous result now? raise IntegralTransformError('Mellin', f, 'could not compute integral') F, cond = F.args[0] if F.has(Integral): raise IntegralTransformError( 'Mellin', f, 'integral in unexpected form') def process_conds(cond): """ Turn ``cond`` into a strip (a, b), and auxiliary conditions. """ from sympy.solvers.inequalities import _solve_inequality a = S.NegativeInfinity b = S.Infinity aux = S.true conds = conjuncts(to_cnf(cond)) t = Dummy('t', real=True) for c in conds: a_ = S.Infinity b_ = S.NegativeInfinity aux_ = [] for d in disjuncts(c): d_ = d.replace( re, lambda x: x.as_real_imag()[0]).subs(re(s), t) if not d.is_Relational or \ d.rel_op in ('==', '!=') \ or d_.has(s) or not d_.has(t): aux_ += [d] continue soln = _solve_inequality(d_, t) if not soln.is_Relational or \ soln.rel_op in ('==', '!='): aux_ += [d] continue if soln.lts == t: b_ = Max(soln.gts, b_) else: a_ = Min(soln.lts, a_) if a_ is not S.Infinity and a_ != b: a = Max(a_, a) elif b_ is not S.NegativeInfinity and b_ != a: b = Min(b_, b) else: aux = And(aux, Or(*aux_)) return a, b, aux conds = [process_conds(c) for c in disjuncts(cond)] conds = [x for x in conds if x[2] != False] conds.sort(key=lambda x: (x[0] - x[1], count_ops(x[2]))) if not conds: raise IntegralTransformError('Mellin', f, 'no convergence found') a, b, aux = conds[0] return _simplify(F.subs(s, s_), simplify), (a, b), aux class MellinTransform(IntegralTransform): """ Class representing unevaluated Mellin transforms. For usage of this class, see the :class:`IntegralTransform` docstring. For how to compute Mellin transforms, see the :func:`mellin_transform` docstring. """ _name = 'Mellin' def _compute_transform(self, f, x, s, **hints): return _mellin_transform(f, x, s, **hints) def _as_integral(self, f, x, s): return Integral(f*x**(s - 1), (x, S.Zero, S.Infinity)) def _collapse_extra(self, extra): a = [] b = [] cond = [] for (sa, sb), c in extra: a += [sa] b += [sb] cond += [c] res = (Max(*a), Min(*b)), And(*cond) if (res[0][0] >= res[0][1]) == True or res[1] == False: raise IntegralTransformError( 'Mellin', None, 'no combined convergence.') return res def mellin_transform(f, x, s, **hints): r""" Compute the Mellin transform `F(s)` of `f(x)`, .. math :: F(s) = \int_0^\infty x^{s-1} f(x) \mathrm{d}x. For all "sensible" functions, this converges absolutely in a strip `a < \operatorname{Re}(s) < b`. Explanation =========== The Mellin transform is related via change of variables to the Fourier transform, and also to the (bilateral) Laplace transform. This function returns ``(F, (a, b), cond)`` where ``F`` is the Mellin transform of ``f``, ``(a, b)`` is the fundamental strip (as above), and ``cond`` are auxiliary convergence conditions. If the integral cannot be computed in closed form, this function returns an unevaluated :class:`MellinTransform` object. For a description of possible hints, refer to the docstring of :func:`sympy.integrals.transforms.IntegralTransform.doit`. If ``noconds=False``, then only `F` will be returned (i.e. not ``cond``, and also not the strip ``(a, b)``). Examples ======== >>> from sympy import mellin_transform, exp >>> from sympy.abc import x, s >>> mellin_transform(exp(-x), x, s) (gamma(s), (0, oo), True) See Also ======== inverse_mellin_transform, laplace_transform, fourier_transform hankel_transform, inverse_hankel_transform """ return MellinTransform(f, x, s).doit(**hints) def _rewrite_sin(m_n, s, a, b): """ Re-write the sine function ``sin(m*s + n)`` as gamma functions, compatible with the strip (a, b). Return ``(gamma1, gamma2, fac)`` so that ``f == fac/(gamma1 * gamma2)``. Examples ======== >>> from sympy.integrals.transforms import _rewrite_sin >>> from sympy import pi, S >>> from sympy.abc import s >>> _rewrite_sin((pi, 0), s, 0, 1) (gamma(s), gamma(1 - s), pi) >>> _rewrite_sin((pi, 0), s, 1, 0) (gamma(s - 1), gamma(2 - s), -pi) >>> _rewrite_sin((pi, 0), s, -1, 0) (gamma(s + 1), gamma(-s), -pi) >>> _rewrite_sin((pi, pi/2), s, S(1)/2, S(3)/2) (gamma(s - 1/2), gamma(3/2 - s), -pi) >>> _rewrite_sin((pi, pi), s, 0, 1) (gamma(s), gamma(1 - s), -pi) >>> _rewrite_sin((2*pi, 0), s, 0, S(1)/2) (gamma(2*s), gamma(1 - 2*s), pi) >>> _rewrite_sin((2*pi, 0), s, S(1)/2, 1) (gamma(2*s - 1), gamma(2 - 2*s), -pi) """ # (This is a separate function because it is moderately complicated, # and I want to doctest it.) # We want to use pi/sin(pi*x) = gamma(x)*gamma(1-x). # But there is one comlication: the gamma functions determine the # inegration contour in the definition of the G-function. Usually # it would not matter if this is slightly shifted, unless this way # we create an undefined function! # So we try to write this in such a way that the gammas are # eminently on the right side of the strip. m, n = m_n m = expand_mul(m/pi) n = expand_mul(n/pi) r = ceiling(-m*a - n.as_real_imag()[0]) # Don't use re(n), does not expand return gamma(m*s + n + r), gamma(1 - n - r - m*s), (-1)**r*pi class MellinTransformStripError(ValueError): """ Exception raised by _rewrite_gamma. Mainly for internal use. """ pass def _rewrite_gamma(f, s, a, b): """ Try to rewrite the product f(s) as a product of gamma functions, so that the inverse Mellin transform of f can be expressed as a meijer G function. Explanation =========== Return (an, ap), (bm, bq), arg, exp, fac such that G((an, ap), (bm, bq), arg/z**exp)*fac is the inverse Mellin transform of f(s). Raises IntegralTransformError or MellinTransformStripError on failure. It is asserted that f has no poles in the fundamental strip designated by (a, b). One of a and b is allowed to be None. The fundamental strip is important, because it determines the inversion contour. This function can handle exponentials, linear factors, trigonometric functions. This is a helper function for inverse_mellin_transform that will not attempt any transformations on f. Examples ======== >>> from sympy.integrals.transforms import _rewrite_gamma >>> from sympy.abc import s >>> from sympy import oo >>> _rewrite_gamma(s*(s+3)*(s-1), s, -oo, oo) (([], [-3, 0, 1]), ([-2, 1, 2], []), 1, 1, -1) >>> _rewrite_gamma((s-1)**2, s, -oo, oo) (([], [1, 1]), ([2, 2], []), 1, 1, 1) Importance of the fundamental strip: >>> _rewrite_gamma(1/s, s, 0, oo) (([1], []), ([], [0]), 1, 1, 1) >>> _rewrite_gamma(1/s, s, None, oo) (([1], []), ([], [0]), 1, 1, 1) >>> _rewrite_gamma(1/s, s, 0, None) (([1], []), ([], [0]), 1, 1, 1) >>> _rewrite_gamma(1/s, s, -oo, 0) (([], [1]), ([0], []), 1, 1, -1) >>> _rewrite_gamma(1/s, s, None, 0) (([], [1]), ([0], []), 1, 1, -1) >>> _rewrite_gamma(1/s, s, -oo, None) (([], [1]), ([0], []), 1, 1, -1) >>> _rewrite_gamma(2**(-s+3), s, -oo, oo) (([], []), ([], []), 1/2, 1, 8) """ # Our strategy will be as follows: # 1) Guess a constant c such that the inversion integral should be # performed wrt s'=c*s (instead of plain s). Write s for s'. # 2) Process all factors, rewrite them independently as gamma functions in # argument s, or exponentials of s. # 3) Try to transform all gamma functions s.t. they have argument # a+s or a-s. # 4) Check that the resulting G function parameters are valid. # 5) Combine all the exponentials. a_, b_ = S([a, b]) def left(c, is_numer): """ Decide whether pole at c lies to the left of the fundamental strip. """ # heuristically, this is the best chance for us to solve the inequalities c = expand(re(c)) if a_ is None and b_ is S.Infinity: return True if a_ is None: return c < b_ if b_ is None: return c <= a_ if (c >= b_) == True: return False if (c <= a_) == True: return True if is_numer: return None if a_.free_symbols or b_.free_symbols or c.free_symbols: return None # XXX #raise IntegralTransformError('Inverse Mellin', f, # 'Could not determine position of singularity %s' # ' relative to fundamental strip' % c) raise MellinTransformStripError('Pole inside critical strip?') # 1) s_multipliers = [] for g in f.atoms(gamma): if not g.has(s): continue arg = g.args[0] if arg.is_Add: arg = arg.as_independent(s)[1] coeff, _ = arg.as_coeff_mul(s) s_multipliers += [coeff] for g in f.atoms(sin, cos, tan, cot): if not g.has(s): continue arg = g.args[0] if arg.is_Add: arg = arg.as_independent(s)[1] coeff, _ = arg.as_coeff_mul(s) s_multipliers += [coeff/pi] s_multipliers = [Abs(x) if x.is_extended_real else x for x in s_multipliers] common_coefficient = S.One for x in s_multipliers: if not x.is_Rational: common_coefficient = x break s_multipliers = [x/common_coefficient for x in s_multipliers] if not (all(x.is_Rational for x in s_multipliers) and common_coefficient.is_extended_real): raise IntegralTransformError("Gamma", None, "Nonrational multiplier") s_multiplier = common_coefficient/reduce(ilcm, [S(x.q) for x in s_multipliers], S.One) if s_multiplier == common_coefficient: if len(s_multipliers) == 0: s_multiplier = common_coefficient else: s_multiplier = common_coefficient \ *reduce(igcd, [S(x.p) for x in s_multipliers]) f = f.subs(s, s/s_multiplier) fac = S.One/s_multiplier exponent = S.One/s_multiplier if a_ is not None: a_ *= s_multiplier if b_ is not None: b_ *= s_multiplier # 2) numer, denom = f.as_numer_denom() numer = Mul.make_args(numer) denom = Mul.make_args(denom) args = list(zip(numer, repeat(True))) + list(zip(denom, repeat(False))) facs = [] dfacs = [] # *_gammas will contain pairs (a, c) representing Gamma(a*s + c) numer_gammas = [] denom_gammas = [] # exponentials will contain bases for exponentials of s exponentials = [] def exception(fact): return IntegralTransformError("Inverse Mellin", f, "Unrecognised form '%s'." % fact) while args: fact, is_numer = args.pop() if is_numer: ugammas, lgammas = numer_gammas, denom_gammas ufacs = facs else: ugammas, lgammas = denom_gammas, numer_gammas ufacs = dfacs def linear_arg(arg): """ Test if arg is of form a*s+b, raise exception if not. """ if not arg.is_polynomial(s): raise exception(fact) p = Poly(arg, s) if p.degree() != 1: raise exception(fact) return p.all_coeffs() # constants if not fact.has(s): ufacs += [fact] # exponentials elif fact.is_Pow or isinstance(fact, exp): if fact.is_Pow: base = fact.base exp_ = fact.exp else: base = exp_polar(1) exp_ = fact.exp if exp_.is_Integer: cond = is_numer if exp_ < 0: cond = not cond args += [(base, cond)]*Abs(exp_) continue elif not base.has(s): a, b = linear_arg(exp_) if not is_numer: base = 1/base exponentials += [base**a] facs += [base**b] else: raise exception(fact) # linear factors elif fact.is_polynomial(s): p = Poly(fact, s) if p.degree() != 1: # We completely factor the poly. For this we need the roots. # Now roots() only works in some cases (low degree), and CRootOf # only works without parameters. So try both... coeff = p.LT()[1] rs = roots(p, s) if len(rs) != p.degree(): rs = CRootOf.all_roots(p) ufacs += [coeff] args += [(s - c, is_numer) for c in rs] continue a, c = p.all_coeffs() ufacs += [a] c /= -a # Now need to convert s - c if left(c, is_numer): ugammas += [(S.One, -c + 1)] lgammas += [(S.One, -c)] else: ufacs += [-1] ugammas += [(S.NegativeOne, c + 1)] lgammas += [(S.NegativeOne, c)] elif isinstance(fact, gamma): a, b = linear_arg(fact.args[0]) if is_numer: if (a > 0 and (left(-b/a, is_numer) == False)) or \ (a < 0 and (left(-b/a, is_numer) == True)): raise NotImplementedError( 'Gammas partially over the strip.') ugammas += [(a, b)] elif isinstance(fact, sin): # We try to re-write all trigs as gammas. This is not in # general the best strategy, since sometimes this is impossible, # but rewriting as exponentials would work. However trig functions # in inverse mellin transforms usually all come from simplifying # gamma terms, so this should work. a = fact.args[0] if is_numer: # No problem with the poles. gamma1, gamma2, fac_ = gamma(a/pi), gamma(1 - a/pi), pi else: gamma1, gamma2, fac_ = _rewrite_sin(linear_arg(a), s, a_, b_) args += [(gamma1, not is_numer), (gamma2, not is_numer)] ufacs += [fac_] elif isinstance(fact, tan): a = fact.args[0] args += [(sin(a, evaluate=False), is_numer), (sin(pi/2 - a, evaluate=False), not is_numer)] elif isinstance(fact, cos): a = fact.args[0] args += [(sin(pi/2 - a, evaluate=False), is_numer)] elif isinstance(fact, cot): a = fact.args[0] args += [(sin(pi/2 - a, evaluate=False), is_numer), (sin(a, evaluate=False), not is_numer)] else: raise exception(fact) fac *= Mul(*facs)/Mul(*dfacs) # 3) an, ap, bm, bq = [], [], [], [] for gammas, plus, minus, is_numer in [(numer_gammas, an, bm, True), (denom_gammas, bq, ap, False)]: while gammas: a, c = gammas.pop() if a != -1 and a != +1: # We use the gamma function multiplication theorem. p = Abs(S(a)) newa = a/p newc = c/p if not a.is_Integer: raise TypeError("a is not an integer") for k in range(p): gammas += [(newa, newc + k/p)] if is_numer: fac *= (2*pi)**((1 - p)/2) * p**(c - S.Half) exponentials += [p**a] else: fac /= (2*pi)**((1 - p)/2) * p**(c - S.Half) exponentials += [p**(-a)] continue if a == +1: plus.append(1 - c) else: minus.append(c) # 4) # TODO # 5) arg = Mul(*exponentials) # for testability, sort the arguments an.sort(key=default_sort_key) ap.sort(key=default_sort_key) bm.sort(key=default_sort_key) bq.sort(key=default_sort_key) return (an, ap), (bm, bq), arg, exponent, fac @_noconds_(True) def _inverse_mellin_transform(F, s, x_, strip, as_meijerg=False): """ A helper for the real inverse_mellin_transform function, this one here assumes x to be real and positive. """ x = _dummy('t', 'inverse-mellin-transform', F, positive=True) # Actually, we won't try integration at all. Instead we use the definition # of the Meijer G function as a fairly general inverse mellin transform. F = F.rewrite(gamma) for g in [factor(F), expand_mul(F), expand(F)]: if g.is_Add: # do all terms separately ress = [_inverse_mellin_transform(G, s, x, strip, as_meijerg, noconds=False) for G in g.args] conds = [p[1] for p in ress] ress = [p[0] for p in ress] res = Add(*ress) if not as_meijerg: res = factor(res, gens=res.atoms(Heaviside)) return res.subs(x, x_), And(*conds) try: a, b, C, e, fac = _rewrite_gamma(g, s, strip[0], strip[1]) except IntegralTransformError: continue try: G = meijerg(a, b, C/x**e) except ValueError: continue if as_meijerg: h = G else: try: from sympy.simplify import hyperexpand h = hyperexpand(G) except NotImplementedError: raise IntegralTransformError( 'Inverse Mellin', F, 'Could not calculate integral') if h.is_Piecewise and len(h.args) == 3: # XXX we break modularity here! h = Heaviside(x - Abs(C))*h.args[0].args[0] \ + Heaviside(Abs(C) - x)*h.args[1].args[0] # We must ensure that the integral along the line we want converges, # and return that value. # See [L], 5.2 cond = [Abs(arg(G.argument)) < G.delta*pi] # Note: we allow ">=" here, this corresponds to convergence if we let # limits go to oo symmetrically. ">" corresponds to absolute convergence. cond += [And(Or(len(G.ap) != len(G.bq), 0 >= re(G.nu) + 1), Abs(arg(G.argument)) == G.delta*pi)] cond = Or(*cond) if cond == False: raise IntegralTransformError( 'Inverse Mellin', F, 'does not converge') return (h*fac).subs(x, x_), cond raise IntegralTransformError('Inverse Mellin', F, '') _allowed = None class InverseMellinTransform(IntegralTransform): """ Class representing unevaluated inverse Mellin transforms. For usage of this class, see the :class:`IntegralTransform` docstring. For how to compute inverse Mellin transforms, see the :func:`inverse_mellin_transform` docstring. """ _name = 'Inverse Mellin' _none_sentinel = Dummy('None') _c = Dummy('c') def __new__(cls, F, s, x, a, b, **opts): if a is None: a = InverseMellinTransform._none_sentinel if b is None: b = InverseMellinTransform._none_sentinel return IntegralTransform.__new__(cls, F, s, x, a, b, **opts) @property def fundamental_strip(self): a, b = self.args[3], self.args[4] if a is InverseMellinTransform._none_sentinel: a = None if b is InverseMellinTransform._none_sentinel: b = None return a, b def _compute_transform(self, F, s, x, **hints): # IntegralTransform's doit will cause this hint to exist, but # InverseMellinTransform should ignore it hints.pop('simplify', True) global _allowed if _allowed is None: _allowed = { exp, gamma, sin, cos, tan, cot, cosh, sinh, tanh, coth, factorial, rf} for f in postorder_traversal(F): if f.is_Function and f.has(s) and f.func not in _allowed: raise IntegralTransformError('Inverse Mellin', F, 'Component %s not recognised.' % f) strip = self.fundamental_strip return _inverse_mellin_transform(F, s, x, strip, **hints) def _as_integral(self, F, s, x): c = self.__class__._c return Integral(F*x**(-s), (s, c - S.ImaginaryUnit*S.Infinity, c + S.ImaginaryUnit*S.Infinity))/(2*S.Pi*S.ImaginaryUnit) def inverse_mellin_transform(F, s, x, strip, **hints): r""" Compute the inverse Mellin transform of `F(s)` over the fundamental strip given by ``strip=(a, b)``. Explanation =========== This can be defined as .. math:: f(x) = \frac{1}{2\pi i} \int_{c - i\infty}^{c + i\infty} x^{-s} F(s) \mathrm{d}s, for any `c` in the fundamental strip. Under certain regularity conditions on `F` and/or `f`, this recovers `f` from its Mellin transform `F` (and vice versa), for positive real `x`. One of `a` or `b` may be passed as ``None``; a suitable `c` will be inferred. If the integral cannot be computed in closed form, this function returns an unevaluated :class:`InverseMellinTransform` object. Note that this function will assume x to be positive and real, regardless of the SymPy assumptions! For a description of possible hints, refer to the docstring of :func:`sympy.integrals.transforms.IntegralTransform.doit`. Examples ======== >>> from sympy import inverse_mellin_transform, oo, gamma >>> from sympy.abc import x, s >>> inverse_mellin_transform(gamma(s), s, x, (0, oo)) exp(-x) The fundamental strip matters: >>> f = 1/(s**2 - 1) >>> inverse_mellin_transform(f, s, x, (-oo, -1)) x*(1 - 1/x**2)*Heaviside(x - 1)/2 >>> inverse_mellin_transform(f, s, x, (-1, 1)) -x*Heaviside(1 - x)/2 - Heaviside(x - 1)/(2*x) >>> inverse_mellin_transform(f, s, x, (1, oo)) (1/2 - x**2/2)*Heaviside(1 - x)/x See Also ======== mellin_transform hankel_transform, inverse_hankel_transform """ return InverseMellinTransform(F, s, x, strip[0], strip[1]).doit(**hints) ########################################################################## # Laplace Transform ########################################################################## def _simplifyconds(expr, s, a): r""" Naively simplify some conditions occurring in ``expr``, given that `\operatorname{Re}(s) > a`. Examples ======== >>> from sympy.integrals.transforms import _simplifyconds as simp >>> from sympy.abc import x >>> from sympy import sympify as S >>> simp(abs(x**2) < 1, x, 1) False >>> simp(abs(x**2) < 1, x, 2) False >>> simp(abs(x**2) < 1, x, 0) Abs(x**2) < 1 >>> simp(abs(1/x**2) < 1, x, 1) True >>> simp(S(1) < abs(x), x, 1) True >>> simp(S(1) < abs(1/x), x, 1) False >>> from sympy import Ne >>> simp(Ne(1, x**3), x, 1) True >>> simp(Ne(1, x**3), x, 2) True >>> simp(Ne(1, x**3), x, 0) Ne(1, x**3) """ def power(ex): if ex == s: return 1 if ex.is_Pow and ex.base == s: return ex.exp return None def bigger(ex1, ex2): """ Return True only if |ex1| > |ex2|, False only if |ex1| < |ex2|. Else return None. """ if ex1.has(s) and ex2.has(s): return None if isinstance(ex1, Abs): ex1 = ex1.args[0] if isinstance(ex2, Abs): ex2 = ex2.args[0] if ex1.has(s): return bigger(1/ex2, 1/ex1) n = power(ex2) if n is None: return None try: if n > 0 and (Abs(ex1) <= Abs(a)**n) == True: return False if n < 0 and (Abs(ex1) >= Abs(a)**n) == True: return True except TypeError: pass def replie(x, y): """ simplify x < y """ if not (x.is_positive or isinstance(x, Abs)) \ or not (y.is_positive or isinstance(y, Abs)): return (x < y) r = bigger(x, y) if r is not None: return not r return (x < y) def replue(x, y): b = bigger(x, y) if b in (True, False): return True return Unequality(x, y) def repl(ex, *args): if ex in (True, False): return bool(ex) return ex.replace(*args) from sympy.simplify.radsimp import collect_abs expr = collect_abs(expr) expr = repl(expr, Lt, replie) expr = repl(expr, Gt, lambda x, y: replie(y, x)) expr = repl(expr, Unequality, replue) return S(expr) def expand_dirac_delta(expr): """ Expand an expression involving DiractDelta to get it as a linear combination of DiracDelta functions. """ return _lin_eq2dict(expr, expr.atoms(DiracDelta)) @_noconds def _laplace_transform(f, t, s_, simplify=True): """ The backend function for Laplace transforms. This backend assumes that the frontend has already split sums such that `f` is to an addition anymore. """ s = Dummy('s') a = Wild('a', exclude=[t]) deltazero = [] deltanonzero = [] try: integratable, deltadict = expand_dirac_delta(f) except PolyNonlinearError: raise IntegralTransformError( 'Laplace', f, 'could not expand DiracDelta expressions') for dirac_func, dirac_coeff in deltadict.items(): p = dirac_func.match(DiracDelta(a*t)) if p: deltazero.append(dirac_coeff.subs(t,0)/p[a]) else: if dirac_func.args[0].subs(t,0).is_zero: raise IntegralTransformError('Laplace', f,\ 'not implemented yet.') else: deltanonzero.append(dirac_func*dirac_coeff) F = Add(integrate(exp(-s*t) * Add(integratable, *deltanonzero), (t, S.Zero, S.Infinity)), Add(*deltazero)) if not F.has(Integral): return _simplify(F.subs(s, s_), simplify), S.NegativeInfinity, S.true if not F.is_Piecewise: raise IntegralTransformError( 'Laplace', f, 'could not compute integral') F, cond = F.args[0] if F.has(Integral): raise IntegralTransformError( 'Laplace', f, 'integral in unexpected form') def process_conds(conds): """ Turn ``conds`` into a strip and auxiliary conditions. """ from sympy.solvers.inequalities import _solve_inequality a = S.NegativeInfinity aux = S.true conds = conjuncts(to_cnf(conds)) p, q, w1, w2, w3, w4, w5 = symbols( 'p q w1 w2 w3 w4 w5', cls=Wild, exclude=[s]) patterns = ( p*Abs(arg((s + w3)*q)) < w2, p*Abs(arg((s + w3)*q)) <= w2, Abs(periodic_argument((s + w3)**p*q, w1)) < w2, Abs(periodic_argument((s + w3)**p*q, w1)) <= w2, Abs(periodic_argument((polar_lift(s + w3))**p*q, w1)) < w2, Abs(periodic_argument((polar_lift(s + w3))**p*q, w1)) <= w2) for c in conds: a_ = S.Infinity aux_ = [] for d in disjuncts(c): if d.is_Relational and s in d.rhs.free_symbols: d = d.reversed if d.is_Relational and isinstance(d, (Ge, Gt)): d = d.reversedsign for pat in patterns: m = d.match(pat) if m: break if m: if m[q].is_positive and m[w2]/m[p] == pi/2: d = -re(s + m[w3]) < 0 m = d.match(p - cos(w1*Abs(arg(s*w5))*w2)*Abs(s**w3)**w4 < 0) if not m: m = d.match( cos(p - Abs(periodic_argument(s**w1*w5, q))*w2)*Abs(s**w3)**w4 < 0) if not m: m = d.match( p - cos(Abs(periodic_argument(polar_lift(s)**w1*w5, q))*w2 )*Abs(s**w3)**w4 < 0) if m and all(m[wild].is_positive for wild in [w1, w2, w3, w4, w5]): d = re(s) > m[p] d_ = d.replace( re, lambda x: x.expand().as_real_imag()[0]).subs(re(s), t) if not d.is_Relational or \ d.rel_op in ('==', '!=') \ or d_.has(s) or not d_.has(t): aux_ += [d] continue soln = _solve_inequality(d_, t) if not soln.is_Relational or \ soln.rel_op in ('==', '!='): aux_ += [d] continue if soln.lts == t: raise IntegralTransformError('Laplace', f, 'convergence not in half-plane?') else: a_ = Min(soln.lts, a_) if a_ is not S.Infinity: a = Max(a_, a) else: aux = And(aux, Or(*aux_)) return a, aux.canonical if aux.is_Relational else aux conds = [process_conds(c) for c in disjuncts(cond)] conds2 = [x for x in conds if x[1] != False and x[0] is not S.NegativeInfinity] if not conds2: conds2 = [x for x in conds if x[1] != False] conds = list(ordered(conds2)) def cnt(expr): if expr in (True, False): return 0 return expr.count_ops() conds.sort(key=lambda x: (-x[0], cnt(x[1]))) if not conds: raise IntegralTransformError('Laplace', f, 'no convergence found') a, aux = conds[0] # XXX is [0] always the right one? def sbs(expr): return expr.subs(s, s_) if simplify: F = _simplifyconds(F, s, a) aux = _simplifyconds(aux, s, a) return _simplify(F.subs(s, s_), simplify), sbs(a), _canonical(sbs(aux)) def _laplace_deep_collect(f, t): """ This is an internal helper function that traverses through the epression tree of `f(t)` and collects arguments. The purpose of it is that anything like `f(w*t-1*t-c)` will be written as `f((w-1)*t-c)` such that it can match `f(a*t+b)`. """ func = f.func args = list(f.args) if len(f.args) == 0: return f else: args = [_laplace_deep_collect(arg, t) for arg in args] if func.is_Add: return func(*args).collect(t) else: return func(*args) def _laplace_build_rules(t, s): """ This is an internal helper function that returns the table of Laplace transfrom rules in terms of the time variable `t` and the frequency variable `s`. It is used by `_laplace_apply_rules`. """ a = Wild('a', exclude=[t]) b = Wild('b', exclude=[t]) n = Wild('n', exclude=[t]) tau = Wild('tau', exclude=[t]) omega = Wild('omega', exclude=[t]) dco = lambda f: _laplace_deep_collect(f,t) laplace_transform_rules = [ # ( time domain, # laplace domain, # condition, convergence plane, preparation function ) # # Catch constant (would otherwise be treated by 2.12) (a, a/s, S.true, S.Zero, dco), # DiracDelta rules (DiracDelta(a*t-b), exp(-s*b/a)/Abs(a), Or(And(a>0, b>=0), And(a<0, b<=0)), S.Zero, dco), (DiracDelta(a*t-b), S(0), Or(And(a<0, b>=0), And(a>0, b<=0)), S.Zero, dco), # Rules from http://eqworld.ipmnet.ru/en/auxiliary/inttrans/ # 2.1 (1, 1/s, S.true, S.Zero, dco), # 2.2 expressed in terms of Heaviside (Heaviside(a*t-b), exp(-s*b/a)/s, And(a>0, b>0), S.Zero, dco), (Heaviside(a*t-b), (1-exp(-s*b/a))/s, And(a<0, b<0), S.Zero, dco), (Heaviside(a*t-b), 1/s, And(a>0, b<=0), S.Zero, dco), (Heaviside(a*t-b), 0, And(a<0, b>0), S.Zero, dco), # 2.3 (t, 1/s**2, S.true, S.Zero, dco), # 2.4 (1/(a*t+b), -exp(-b/a*s)*Ei(-b/a*s)/a, a>0, S.Zero, dco), # 2.5 and 2.6 are covered by 2.11 # 2.7 (1/sqrt(a*t+b), sqrt(a*pi/s)*exp(b/a*s)*erfc(sqrt(b/a*s))/a, a>0, S.Zero, dco), # 2.8 (sqrt(t)/(t+b), sqrt(pi/s)-pi*sqrt(b)*exp(b*s)*erfc(sqrt(b*s)), S.true, S.Zero, dco), # 2.9 ((a*t+b)**(-S(3)/2), 2*b**(-S(1)/2)-2*(pi*s/a)**(S(1)/2)*exp(b/a*s)*erfc(sqrt(b/a*s))/a, a>0, S.Zero, dco), # 2.10 (t**(S(1)/2)*(t+a)**(-1), (pi/s)**(S(1)/2)-pi*a**(S(1)/2)*exp(a*s)*erfc(sqrt(a*s)), S.true, S.Zero, dco), # 2.11 (1/(a*sqrt(t) + t**(3/2)), pi*a**(S(1)/2)*exp(a*s)*erfc(sqrt(a*s)), S.true, S.Zero, dco), # 2.12 (t**n, gamma(n+1)/s**(n+1), n>-1, S.Zero, dco), # 2.13 ((a*t+b)**n, lowergamma(n+1, b/a*s)*exp(-b/a*s)/s**(n+1)/a, And(n>-1, a>0), S.Zero, dco), # 2.14 (t**n/(t+a), a**n*gamma(n+1)*lowergamma(-n,a*s), n>-1, S.Zero, dco), # 3.1 (exp(a*t-tau), exp(-tau)/(s-a), S.true, a, dco), # 3.2 (t*exp(a*t-tau), exp(-tau)/(s-a)**2, S.true, a, dco), # 3.3 (t**n*exp(a*t), gamma(n+1)/(s-a)**(n+1), n>-1, a, dco), # 3.4 and 3.5 cannot be covered here because they are # sums and only the individual sum terms will get here. # 3.6 (exp(-a*t**2), sqrt(pi/4/a)*exp(s**2/4/a)*erfc(s/sqrt(4*a)), a>0, S.Zero, dco), # 3.7 (t*exp(-a*t**2), 1/(2*a)-2/sqrt(pi)/(4*a)**(S(3)/2)*s*erfc(s/sqrt(4*a)), S.true, S.Zero, dco), # 3.8 (exp(-a/t), 2*sqrt(a/s)*besselk(1, 2*sqrt(a*s)), a>=0, S.Zero, dco), # 3.9 (sqrt(t)*exp(-a/t), S(1)/2*sqrt(pi/s**3)*(1+2*sqrt(a*s))*exp(-2*sqrt(a*s)), a>=0, S.Zero, dco), # 3.10 (exp(-a/t)/sqrt(t), sqrt(pi/s)*exp(-2*sqrt(a*s)), a>=0, S.Zero, dco), # 3.11 (exp(-a/t)/(t*sqrt(t)), sqrt(pi/a)*exp(-2*sqrt(a*s)), a>0, S.Zero, dco), # 3.12 (t**n*exp(-a/t), 2*(a/s)**((n+1)/2)*besselk(n+1, 2*sqrt(a*s)), a>0, S.Zero, dco), # 3.13 (exp(-2*sqrt(a*t)), s**(-1)-sqrt(pi*a)*s**(-S(3)/2)*exp(a/s)*erfc(sqrt(a/s)), S.true, S.Zero, dco), # 3.14 (exp(-2*sqrt(a*t))/sqrt(t), (pi/s)**(S(1)/2)*exp(a/s)*erfc(sqrt(a/s)), S.true, S.Zero, dco), # 4.1 (sinh(a*t), a/(s**2-a**2), S.true, Abs(a), dco), # 4.2 (sinh(a*t)**2, 2*a**2/(s**3-4*a**2*s**2), S.true, Abs(2*a), dco), # 4.3 (sinh(a*t)/t, log((s+a)/(s-a))/2, S.true, a, dco), # 4.4 (t**n*sinh(a*t), gamma(n+1)/2*((s-a)**(-n-1)-(s+a)**(-n-1)), n>-2, Abs(a), dco), # 4.5 (sinh(2*sqrt(a*t)), sqrt(pi*a)/s/sqrt(s)*exp(a/s), S.true, S.Zero, dco), # 4.6 (sqrt(t)*sinh(2*sqrt(a*t)), pi**(S(1)/2)*s**(-S(5)/2)*(s/2+a)*exp(a/s)*erf(sqrt(a/s))-a**(S(1)/2)*s**(-2), S.true, S.Zero, dco), # 4.7 (sinh(2*sqrt(a*t))/sqrt(t), pi**(S(1)/2)*s**(-S(1)/2)*exp(a/s)*erf(sqrt(a/s)), S.true, S.Zero, dco), # 4.8 (sinh(sqrt(a*t))**2/sqrt(t), pi**(S(1)/2)/2*s**(-S(1)/2)*(exp(a/s)-1), S.true, S.Zero, dco), # 4.9 (cosh(a*t), s/(s**2-a**2), S.true, Abs(a), dco), # 4.10 (cosh(a*t)**2, (s**2-2*a**2)/(s**3-4*a**2*s**2), S.true, Abs(2*a), dco), # 4.11 (t**n*cosh(a*t), gamma(n+1)/2*((s-a)**(-n-1)+(s+a)**(-n-1)), n>-1, Abs(a), dco), # 4.12 (cosh(2*sqrt(a*t)), 1/s+sqrt(pi*a)/s/sqrt(s)*exp(a/s)*erf(sqrt(a/s)), S.true, S.Zero, dco), # 4.13 (sqrt(t)*cosh(2*sqrt(a*t)), pi**(S(1)/2)*s**(-S(5)/2)*(s/2+a)*exp(a/s), S.true, S.Zero, dco), # 4.14 (cosh(2*sqrt(a*t))/sqrt(t), pi**(S(1)/2)*s**(-S(1)/2)*exp(a/s), S.true, S.Zero, dco), # 4.15 (cosh(sqrt(a*t))**2/sqrt(t), pi**(S(1)/2)/2*s**(-S(1)/2)*(exp(a/s)+1), S.true, S.Zero, dco), # 5.1 (log(a*t), -log(s/a+S.EulerGamma)/s, a>0, S.Zero, dco), # 5.2 (log(1+a*t), -exp(s/a)/s*Ei(-s/a), S.true, S.Zero, dco), # 5.3 (log(a*t+b), (log(b)-exp(s/b/a)/s*a*Ei(-s/b))/s*a, a>0, S.Zero, dco), # 5.4 is covered by 5.7 # 5.5 (log(t)/sqrt(t), -sqrt(pi/s)*(log(4*s)+S.EulerGamma), S.true, S.Zero, dco), # 5.6 is covered by 5.7 # 5.7 (t**n*log(t), gamma(n+1)*s**(-n-1)*(digamma(n+1)-log(s)), n>-1, S.Zero, dco), # 5.8 (log(a*t)**2, ((log(s/a)+S.EulerGamma)**2+pi**2/6)/s, a>0, S.Zero, dco), # 5.9 (exp(-a*t)*log(t), -(log(s+a)+S.EulerGamma)/(s+a), S.true, -a, dco), # 6.1 (sin(omega*t), omega/(s**2+omega**2), S.true, S.Zero, dco), # 6.2 (Abs(sin(omega*t)), omega/(s**2+omega**2)*coth(pi*s/2/omega), omega>0, S.Zero, dco), # 6.3 and 6.4 are covered by 1.8 # 6.5 is covered by 1.8 together with 2.5 # 6.6 (sin(omega*t)/t, atan(omega/s), S.true, S.Zero, dco), # 6.7 (sin(omega*t)**2/t, log(1+4*omega**2/s**2)/4, S.true, S.Zero, dco), # 6.8 (sin(omega*t)**2/t**2, omega*atan(2*omega/s)-s*log(1+4*omega**2/s**2)/4, S.true, S.Zero, dco), # 6.9 (sin(2*sqrt(a*t)), sqrt(pi*a)/s/sqrt(s)*exp(-a/s), a>0, S.Zero, dco), # 6.10 (sin(2*sqrt(a*t))/t, pi*erf(sqrt(a/s)), a>0, S.Zero, dco), # 6.11 (cos(omega*t), s/(s**2+omega**2), S.true, S.Zero, dco), # 6.12 (cos(omega*t)**2, (s**2+2*omega**2)/(s**2+4*omega**2)/s, S.true, S.Zero, dco), # 6.13 is covered by 1.9 together with 2.5 # 6.14 and 6.15 cannot be done with this method, the respective sum # parts do not converge. Solve elsewhere if really needed. # 6.16 (sqrt(t)*cos(2*sqrt(a*t)), sqrt(pi)/2*s**(-S(5)/2)*(s-2*a)*exp(-a/s), a>0, S.Zero, dco), # 6.17 (cos(2*sqrt(a*t))/sqrt(t), sqrt(pi/s)*exp(-a/s), a>0, S.Zero, dco), # 6.18 (sin(a*t)*sin(b*t), 2*a*b*s/(s**2+(a+b)**2)/(s**2+(a-b)**2), S.true, S.Zero, dco), # 6.19 (cos(a*t)*sin(b*t), b*(s**2-a**2+b**2)/(s**2+(a+b)**2)/(s**2+(a-b)**2), S.true, S.Zero, dco), # 6.20 (cos(a*t)*cos(b*t), s*(s**2+a**2+b**2)/(s**2+(a+b)**2)/(s**2+(a-b)**2), S.true, S.Zero, dco), # 6.21 (exp(b*t)*sin(a*t), a/((s-b)**2+a**2), S.true, b, dco), # 6.22 (exp(b*t)*cos(a*t), (s-b)/((s-b)**2+a**2), S.true, b, dco), # 7.1 (erf(a*t), exp(s**2/(2*a)**2)*erfc(s/(2*a))/s, a>0, S.Zero, dco), # 7.2 (erf(sqrt(a*t)), sqrt(a)/sqrt(s+a)/s, a>0, S.Zero, dco), # 7.3 (exp(a*t)*erf(sqrt(a*t)), sqrt(a)/sqrt(s)/(s-a), a>0, a, dco), # 7.4 (erf(sqrt(a/t)/2), (1-exp(-sqrt(a*s)))/s, a>0, S.Zero, dco), # 7.5 (erfc(sqrt(a*t)), (sqrt(s+a)-sqrt(a))/sqrt(s+a)/s, a>0, S.Zero, dco), # 7.6 (exp(a*t)*erfc(sqrt(a*t)), 1/(s+sqrt(a*s)), a>0, S.Zero, dco), # 7.7 (erfc(sqrt(a/t)/2), exp(-sqrt(a*s))/s, a>0, S.Zero, dco), # 8.1, 8.2 (besselj(n, a*t), a**n/(sqrt(s**2+a**2)*(s+sqrt(s**2+a**2))**n), And(a>0, n>-1), S.Zero, dco), # 8.3, 8.4 (t**b*besselj(n, a*t), 2**n/sqrt(pi)*gamma(n+S.Half)*a**n*(s**2+a**2)**(-n-S.Half), And(And(a>0, n>-S.Half), Eq(b, n)), S.Zero, dco), # 8.5 (t**b*besselj(n, a*t), 2**(n+1)/sqrt(pi)*gamma(n+S(3)/2)*a**n*s*(s**2+a**2)**(-n-S(3)/2), And(And(a>0, n>-1), Eq(b, n+1)), S.Zero, dco), # 8.6 (besselj(0, 2*sqrt(a*t)), exp(-a/s)/s, a>0, S.Zero, dco), # 8.7, 8.8 (t**(b)*besselj(n, 2*sqrt(a*t)), a**(n/2)*s**(-n-1)*exp(-a/s), And(And(a>0, n>-1), Eq(b, n*S.Half)), S.Zero, dco), # 8.9 (besselj(0, a*sqrt(t**2+b*t)), exp(b*s-b*sqrt(s**2+a**2))/sqrt(s**2+a**2), b>0, S.Zero, dco), # 8.10, 8.11 (besseli(n, a*t), a**n/(sqrt(s**2-a**2)*(s+sqrt(s**2-a**2))**n), And(a>0, n>-1), Abs(a), dco), # 8.12 (t**b*besseli(n, a*t), 2**n/sqrt(pi)*gamma(n+S.Half)*a**n*(s**2-a**2)**(-n-S.Half), And(And(a>0, n>-S.Half), Eq(b, n)), Abs(a), dco), # 8.13 (t**b*besseli(n, a*t), 2**(n+1)/sqrt(pi)*gamma(n+S(3)/2)*a**n*s*(s**2-a**2)**(-n-S(3)/2), And(And(a>0, n>-1), Eq(b, n+1)), Abs(a), dco), # 8.15, 8.16 (t**(b)*besseli(n, 2*sqrt(a*t)), a**(n/2)*s**(-n-1)*exp(a/s), And(And(a>0, n>-1), Eq(b, n*S.Half)), S.Zero, dco), # 8.17 (bessely(0, a*t), -2/pi*asinh(s/a)/sqrt(s**2+a**2), a>0, S.Zero, dco), # 8.18 (besselk(0, a*t), (log(s+sqrt(s**2-a**2)))/(sqrt(s**2-a**2)), a>0, Abs(a), dco) ] return laplace_transform_rules def _laplace_cr(f, a, c, **hints): """ Internal helper function that will return `(f, a, c)` unless `**hints` contains `noconds=True`, in which case it will only return `f`. """ conds = not hints.get('noconds', False) if conds: return f, a, c else: return f def _laplace_rule_timescale(f, t, s, doit=True, **hints): r""" This internal helper function tries to apply the time-scaling rule of the Laplace transform and returns `None` if it cannot do it. Time-scaling means the following: if $F(s)$ is the Laplace transform of, $f(t)$, then, for any $a>0$, the Laplace transform of $f(at)$ will be $\frac1a F(\frac{s}{a})$. This scaling will also affect the transform's convergence plane. """ _simplify = hints.pop('simplify', True) b = Wild('b', exclude=[t]) g = WildFunction('g', nargs=1) k, func = f.as_independent(t, as_Add=False) ma1 = func.match(g) if ma1: arg = ma1[g].args[0].collect(t) ma2 = arg.match(b*t) if ma2 and ma2[b]>0: debug('_laplace_apply_rules match:') debug(' f: %s ( %s, %s )'%(f, ma1, ma2)) debug(' rule: amplitude and time scaling (1.1, 1.2)') if ma2[b]==1: if doit==True and not any(func.has(t) for func in ma1[g].atoms(AppliedUndef)): return k*_laplace_transform(ma1[g].func(t), t, s, simplify=_simplify) else: return k*LaplaceTransform(ma1[g].func(t), t, s, **hints) else: L = _laplace_apply_rules(ma1[g].func(t), t, s/ma2[b], doit=doit, **hints) try: r, p, c = L return (k/ma2[b]*r, p, c) except TypeError: return k/ma2[b]*L return None def _laplace_rule_heaviside(f, t, s, doit=True, **hints): """ This internal helper function tries to transform a product containing the `Heaviside` function and returns `None` if it cannot do it. """ hints.pop('simplify', True) a = Wild('a', exclude=[t]) b = Wild('b', exclude=[t]) y = Wild('y') g = WildFunction('g', nargs=1) k, func = f.as_independent(t, as_Add=False) ma1 = func.match(Heaviside(y)*g) if ma1: ma2 = ma1[y].match(t-a) ma3 = ma1[g].args[0].collect(t).match(t-b) if ma2 and ma2[a]>0 and ma3 and ma2[a]==ma3[b]: debug('_laplace_apply_rules match:') debug(' f: %s ( %s, %s, %s )'%(f, ma1, ma2, ma3)) debug(' rule: time shift (1.3)') L = _laplace_apply_rules(ma1[g].func(t), t, s, doit=doit, **hints) try: r, p, c = L return (k*exp(-ma2[a]*s)*r, p, c) except TypeError: return k*exp(-ma2[a]*s)*L return None def _laplace_rule_exp(f, t, s, doit=True, **hints): """ This internal helper function tries to transform a product containing the `exp` function and returns `None` if it cannot do it. """ hints.pop('simplify', True) a = Wild('a', exclude=[t]) y = Wild('y') z = Wild('z') k, func = f.as_independent(t, as_Add=False) ma1 = func.match(exp(y)*z) if ma1: ma2 = ma1[y].collect(t).match(a*t) if ma2: debug('_laplace_apply_rules match:') debug(' f: %s ( %s, %s )'%(f, ma1, ma2)) debug(' rule: multiply with exp (1.5)') L = _laplace_apply_rules(ma1[z], t, s-ma2[a], doit=doit, **hints) try: r, p, c = L return (r, p+ma2[a], c) except TypeError: return L return None def _laplace_rule_trig(f, t, s, doit=True, **hints): """ This internal helper function tries to transform a product containing a trigonometric function (`sin`, `cos`, `sinh`, `cosh`, ) and returns `None` if it cannot do it. """ _simplify = hints.pop('simplify', True) a = Wild('a', exclude=[t]) y = Wild('y') z = Wild('z') k, func = f.as_independent(t, as_Add=False) # All of the rules have a very similar form: trig(y)*z is matched, and then # two copies of the Laplace transform of z are shifted in the s Domain # and added with a weight; see rules 1.6 to 1.9 in # http://eqworld.ipmnet.ru/en/auxiliary/inttrans/laplace1.pdf # The parameters in the tuples are (fm, nu, s1, s2, sd): # fm: Function to match # nu: Number of the rule, for debug purposes # s1: weight of the sum, 'I' for sin and '1' for all others # s2: sign of the second copy of the Laplace transform of z # sd: shift direction; shift along real or imaginary axis if `1` or `I` trigrules = [(sinh(y), '1.6', 1, -1, 1), (cosh(y), '1.7', 1, 1, 1), (sin(y), '1.8', -I, -1, I), (cos(y), '1.9', 1, 1, I)] for trigrule in trigrules: fm, nu, s1, s2, sd = trigrule ma1 = func.match(fm*z) if ma1: ma2 = ma1[y].collect(t).match(a*t) if ma2: debug('_laplace_apply_rules match:') debug(' f: %s ( %s, %s )'%(f, ma1, ma2)) debug(' rule: multiply with %s (%s)'%(fm.func, nu)) L = _laplace_apply_rules(ma1[z], t, s, doit=doit, **hints) try: r, p, c = L # The convergence plane changes only if the shift has been # done along the real axis: if sd==1: cp_shift = Abs(ma2[a]) else: cp_shift = 0 return ((s1*(r.subs(s, s-sd*ma2[a])+\ s2*r.subs(s, s+sd*ma2[a]))).simplify()/2, p+cp_shift, c) except TypeError: if doit==True and _simplify==True: return (s1*(L.subs(s, s-sd*ma2[a])+\ s2*L.subs(s, s+sd*ma2[a]))).simplify()/2 else: return (s1*(L.subs(s, s-sd*ma2[a])+\ s2*L.subs(s, s+sd*ma2[a])))/2 return None def _laplace_rule_diff(f, t, s, doit=True, **hints): """ This internal helper function tries to transform an expression containing a derivative of an undefined function and returns `None` if it cannot do it. """ hints.pop('simplify', True) a = Wild('a', exclude=[t]) y = Wild('y') n = Wild('n', exclude=[t]) g = WildFunction('g', nargs=1) ma1 = f.match(a*Derivative(g, (t, n))) if ma1 and ma1[g].args[0] == t and ma1[n].is_integer: debug('_laplace_apply_rules match:') debug(' f: %s'%(f,)) debug(' rule: time derivative (1.11, 1.12)') d = [] for k in range(ma1[n]): if k==0: y = ma1[g].func(t).subs(t, 0) else: y = Derivative(ma1[g].func(t), (t, k)).subs(t, 0) d.append(s**(ma1[n]-k-1)*y) r = s**ma1[n]*_laplace_apply_rules(ma1[g].func(t), t, s, doit=doit, **hints) return ma1[a]*(r - Add(*d)) return None def _laplace_apply_rules(f, t, s, doit=True, **hints): """ Helper function for the class LaplaceTransform. This function does a Laplace transform based on rules and, after applying the rules, hands the rest over to `_laplace_transform`, which will attempt to integrate. If it is called with `doit=False`, then it will instead return `LaplaceTransform` objects. """ k, func = f.as_independent(t, as_Add=False) simple_rules = _laplace_build_rules(t, s) for t_dom, s_dom, check, plane, prep in simple_rules: ma = prep(func).match(t_dom) if ma: debug('_laplace_apply_rules match:') debug(' f: %s'%(func,)) debug(' rule: %s o---o %s'%(t_dom, s_dom)) try: debug(' try %s'%(check,)) c = check.xreplace(ma) debug(' check %s -> %s'%(check, c)) if c==True: return _laplace_cr(k*s_dom.xreplace(ma), plane.xreplace(ma), S.true, **hints) except Exception: debug('_laplace_apply_rules did not match.') if f.has(DiracDelta): return None prog_rules = [_laplace_rule_timescale, _laplace_rule_heaviside, _laplace_rule_exp, _laplace_rule_trig, _laplace_rule_diff] for p_rule in prog_rules: LT = p_rule(f, t, s, doit=doit, **hints) if LT is not None: return LT return None class LaplaceTransform(IntegralTransform): """ Class representing unevaluated Laplace transforms. For usage of this class, see the :class:`IntegralTransform` docstring. For how to compute Laplace transforms, see the :func:`laplace_transform` docstring. """ _name = 'Laplace' def _compute_transform(self, f, t, s, **hints): LT = _laplace_apply_rules(f, t, s, **hints) if LT is None: _simplify = hints.pop('simplify', True) debug('_laplace_apply_rules could not match function %s'%(f,)) debug(' hints: %s'%(hints,)) return _laplace_transform(f, t, s, simplify=_simplify, **hints) else: return LT def _as_integral(self, f, t, s): return Integral(f*exp(-s*t), (t, S.Zero, S.Infinity)) def _collapse_extra(self, extra): conds = [] planes = [] for plane, cond in extra: conds.append(cond) planes.append(plane) cond = And(*conds) plane = Max(*planes) if cond == False: raise IntegralTransformError( 'Laplace', None, 'No combined convergence.') return plane, cond def _try_directly(self, **hints): fn = self.function debug('----> _try_directly: %s'%(fn, )) t_ = self.function_variable s_ = self.transform_variable LT = None if not fn.is_Add: fn = expand_mul(fn) try: LT = self._compute_transform(fn, t_, s_, **hints) except IntegralTransformError: LT = None return fn, LT def laplace_transform(f, t, s, legacy_matrix=True, **hints): r""" Compute the Laplace Transform `F(s)` of `f(t)`, .. math :: F(s) = \int_{0^{-}}^\infty e^{-st} f(t) \mathrm{d}t. Explanation =========== For all sensible functions, this converges absolutely in a half-plane .. math :: a < \operatorname{Re}(s) This function returns ``(F, a, cond)`` where ``F`` is the Laplace transform of ``f``, `a` is the half-plane of convergence, and `cond` are auxiliary convergence conditions. The implementation is rule-based, and if you are interested in which rules are applied, and whether integration is attemped, you can switch debug information on by setting ``sympy.SYMPY_DEBUG=True``. The lower bound is `0-`, meaning that this bound should be approached from the lower side. This is only necessary if distributions are involved. At present, it is only done if `f(t)` contains ``DiracDelta``, in which case the Laplace transform is computed implicitly as .. math :: F(s) = \lim_{\tau\to 0^{-}} \int_{\tau}^\infty e^{-st} f(t) \mathrm{d}t by applying rules. If the integral cannot be fully computed in closed form, this function returns an unevaluated :class:`LaplaceTransform` object. For a description of possible hints, refer to the docstring of :func:`sympy.integrals.transforms.IntegralTransform.doit`. If ``noconds=True``, only `F` will be returned (i.e. not ``cond``, and also not the plane ``a``). .. deprecated:: 1.9 Legacy behavior for matrices where ``laplace_transform`` with ``noconds=False`` (the default) returns a Matrix whose elements are tuples. The behavior of ``laplace_transform`` for matrices will change in a future release of SymPy to return a tuple of the transformed Matrix and the convergence conditions for the matrix as a whole. Use ``legacy_matrix=False`` to enable the new behavior. Examples ======== >>> from sympy import DiracDelta, exp, laplace_transform >>> from sympy.abc import t, s, a >>> laplace_transform(t**4, t, s) (24/s**5, 0, True) >>> laplace_transform(t**a, t, s) (gamma(a + 1)/(s*s**a), 0, re(a) > -1) >>> laplace_transform(DiracDelta(t)-a*exp(-a*t),t,s) (s/(a + s), Max(0, -a), True) See Also ======== inverse_laplace_transform, mellin_transform, fourier_transform hankel_transform, inverse_hankel_transform """ debug('\n***** laplace_transform(%s, %s, %s)'%(f, t, s)) if isinstance(f, MatrixBase) and hasattr(f, 'applyfunc'): conds = not hints.get('noconds', False) if conds and legacy_matrix: sympy_deprecation_warning( """ Calling laplace_transform() on a Matrix with noconds=False (the default) is deprecated. Either noconds=True or use legacy_matrix=False to get the new behavior. """, deprecated_since_version="1.9", active_deprecations_target="deprecated-laplace-transform-matrix", ) # Temporarily disable the deprecation warning for non-Expr objects # in Matrix with ignore_warnings(SymPyDeprecationWarning): return f.applyfunc(lambda fij: laplace_transform(fij, t, s, **hints)) else: elements_trans = [laplace_transform(fij, t, s, **hints) for fij in f] if conds: elements, avals, conditions = zip(*elements_trans) f_laplace = type(f)(*f.shape, elements) return f_laplace, Max(*avals), And(*conditions) else: return type(f)(*f.shape, elements_trans) return LaplaceTransform(f, t, s).doit(**hints) @_noconds_(True) def _inverse_laplace_transform(F, s, t_, plane, simplify=True): """ The backend function for inverse Laplace transforms. """ from sympy.integrals.meijerint import meijerint_inversion, _get_coeff_exp # There are two strategies we can try: # 1) Use inverse mellin transforms - related by a simple change of variables. # 2) Use the inversion integral. t = Dummy('t', real=True) def pw_simp(*args): """ Simplify a piecewise expression from hyperexpand. """ # XXX we break modularity here! if len(args) != 3: return Piecewise(*args) arg = args[2].args[0].argument coeff, exponent = _get_coeff_exp(arg, t) e1 = args[0].args[0] e2 = args[1].args[0] return Heaviside(1/Abs(coeff) - t**exponent)*e1 \ + Heaviside(t**exponent - 1/Abs(coeff))*e2 if F.is_rational_function(s): F = F.apart(s) if F.is_Add: f = Add(*[_inverse_laplace_transform(X, s, t, plane, simplify)\ for X in F.args]) return _simplify(f.subs(t, t_), simplify), True try: f, cond = inverse_mellin_transform(F, s, exp(-t), (None, S.Infinity), needeval=True, noconds=False) except IntegralTransformError: f = None if f is None: f = meijerint_inversion(F, s, t) if f is None: raise IntegralTransformError('Inverse Laplace', f, '') if f.is_Piecewise: f, cond = f.args[0] if f.has(Integral): raise IntegralTransformError('Inverse Laplace', f, 'inversion integral of unrecognised form.') else: cond = S.true f = f.replace(Piecewise, pw_simp) if f.is_Piecewise: # many of the functions called below can't work with piecewise # (b/c it has a bool in args) return f.subs(t, t_), cond u = Dummy('u') def simp_heaviside(arg, H0=S.Half): a = arg.subs(exp(-t), u) if a.has(t): return Heaviside(arg, H0) from sympy.solvers.inequalities import _solve_inequality rel = _solve_inequality(a > 0, u) if rel.lts == u: k = log(rel.gts) return Heaviside(t + k, H0) else: k = log(rel.lts) return Heaviside(-(t + k), H0) f = f.replace(Heaviside, simp_heaviside) def simp_exp(arg): return expand_complex(exp(arg)) f = f.replace(exp, simp_exp) # TODO it would be nice to fix cosh and sinh ... simplify messes these # exponentials up return _simplify(f.subs(t, t_), simplify), cond class InverseLaplaceTransform(IntegralTransform): """ Class representing unevaluated inverse Laplace transforms. For usage of this class, see the :class:`IntegralTransform` docstring. For how to compute inverse Laplace transforms, see the :func:`inverse_laplace_transform` docstring. """ _name = 'Inverse Laplace' _none_sentinel = Dummy('None') _c = Dummy('c') def __new__(cls, F, s, x, plane, **opts): if plane is None: plane = InverseLaplaceTransform._none_sentinel return IntegralTransform.__new__(cls, F, s, x, plane, **opts) @property def fundamental_plane(self): plane = self.args[3] if plane is InverseLaplaceTransform._none_sentinel: plane = None return plane def _compute_transform(self, F, s, t, **hints): return _inverse_laplace_transform(F, s, t, self.fundamental_plane, **hints) def _as_integral(self, F, s, t): c = self.__class__._c return Integral(exp(s*t)*F, (s, c - S.ImaginaryUnit*S.Infinity, c + S.ImaginaryUnit*S.Infinity))/(2*S.Pi*S.ImaginaryUnit) def inverse_laplace_transform(F, s, t, plane=None, **hints): r""" Compute the inverse Laplace transform of `F(s)`, defined as .. math :: f(t) = \frac{1}{2\pi i} \int_{c-i\infty}^{c+i\infty} e^{st} F(s) \mathrm{d}s, for `c` so large that `F(s)` has no singularites in the half-plane `\operatorname{Re}(s) > c-\epsilon`. Explanation =========== The plane can be specified by argument ``plane``, but will be inferred if passed as None. Under certain regularity conditions, this recovers `f(t)` from its Laplace Transform `F(s)`, for non-negative `t`, and vice versa. If the integral cannot be computed in closed form, this function returns an unevaluated :class:`InverseLaplaceTransform` object. Note that this function will always assume `t` to be real, regardless of the SymPy assumption on `t`. For a description of possible hints, refer to the docstring of :func:`sympy.integrals.transforms.IntegralTransform.doit`. Examples ======== >>> from sympy import inverse_laplace_transform, exp, Symbol >>> from sympy.abc import s, t >>> a = Symbol('a', positive=True) >>> inverse_laplace_transform(exp(-a*s)/s, s, t) Heaviside(-a + t) See Also ======== laplace_transform, _fast_inverse_laplace hankel_transform, inverse_hankel_transform """ if isinstance(F, MatrixBase) and hasattr(F, 'applyfunc'): return F.applyfunc(lambda Fij: inverse_laplace_transform(Fij, s, t, plane, **hints)) return InverseLaplaceTransform(F, s, t, plane).doit(**hints) def _fast_inverse_laplace(e, s, t): """Fast inverse Laplace transform of rational function including RootSum""" a, b, n = symbols('a, b, n', cls=Wild, exclude=[s]) def _ilt(e): if not e.has(s): return e elif e.is_Add: return _ilt_add(e) elif e.is_Mul: return _ilt_mul(e) elif e.is_Pow: return _ilt_pow(e) elif isinstance(e, RootSum): return _ilt_rootsum(e) else: raise NotImplementedError def _ilt_add(e): return e.func(*map(_ilt, e.args)) def _ilt_mul(e): coeff, expr = e.as_independent(s) if expr.is_Mul: raise NotImplementedError return coeff * _ilt(expr) def _ilt_pow(e): match = e.match((a*s + b)**n) if match is not None: nm, am, bm = match[n], match[a], match[b] if nm.is_Integer and nm < 0: return t**(-nm-1)*exp(-(bm/am)*t)/(am**-nm*gamma(-nm)) if nm == 1: return exp(-(bm/am)*t) / am raise NotImplementedError def _ilt_rootsum(e): expr = e.fun.expr [variable] = e.fun.variables return RootSum(e.poly, Lambda(variable, together(_ilt(expr)))) return _ilt(e) ########################################################################## # Fourier Transform ########################################################################## @_noconds_(True) def _fourier_transform(f, x, k, a, b, name, simplify=True): r""" Compute a general Fourier-type transform .. math:: F(k) = a \int_{-\infty}^{\infty} e^{bixk} f(x)\, dx. For suitable choice of *a* and *b*, this reduces to the standard Fourier and inverse Fourier transforms. """ F = integrate(a*f*exp(b*S.ImaginaryUnit*x*k), (x, S.NegativeInfinity, S.Infinity)) if not F.has(Integral): return _simplify(F, simplify), S.true integral_f = integrate(f, (x, S.NegativeInfinity, S.Infinity)) if integral_f in (S.NegativeInfinity, S.Infinity, S.NaN) or integral_f.has(Integral): raise IntegralTransformError(name, f, 'function not integrable on real axis') if not F.is_Piecewise: raise IntegralTransformError(name, f, 'could not compute integral') F, cond = F.args[0] if F.has(Integral): raise IntegralTransformError(name, f, 'integral in unexpected form') return _simplify(F, simplify), cond class FourierTypeTransform(IntegralTransform): """ Base class for Fourier transforms.""" def a(self): raise NotImplementedError( "Class %s must implement a(self) but does not" % self.__class__) def b(self): raise NotImplementedError( "Class %s must implement b(self) but does not" % self.__class__) def _compute_transform(self, f, x, k, **hints): return _fourier_transform(f, x, k, self.a(), self.b(), self.__class__._name, **hints) def _as_integral(self, f, x, k): a = self.a() b = self.b() return Integral(a*f*exp(b*S.ImaginaryUnit*x*k), (x, S.NegativeInfinity, S.Infinity)) class FourierTransform(FourierTypeTransform): """ Class representing unevaluated Fourier transforms. For usage of this class, see the :class:`IntegralTransform` docstring. For how to compute Fourier transforms, see the :func:`fourier_transform` docstring. """ _name = 'Fourier' def a(self): return 1 def b(self): return -2*S.Pi def fourier_transform(f, x, k, **hints): r""" Compute the unitary, ordinary-frequency Fourier transform of ``f``, defined as .. math:: F(k) = \int_{-\infty}^\infty f(x) e^{-2\pi i x k} \mathrm{d} x. Explanation =========== If the transform cannot be computed in closed form, this function returns an unevaluated :class:`FourierTransform` object. For other Fourier transform conventions, see the function :func:`sympy.integrals.transforms._fourier_transform`. For a description of possible hints, refer to the docstring of :func:`sympy.integrals.transforms.IntegralTransform.doit`. Note that for this transform, by default ``noconds=True``. Examples ======== >>> from sympy import fourier_transform, exp >>> from sympy.abc import x, k >>> fourier_transform(exp(-x**2), x, k) sqrt(pi)*exp(-pi**2*k**2) >>> fourier_transform(exp(-x**2), x, k, noconds=False) (sqrt(pi)*exp(-pi**2*k**2), True) See Also ======== inverse_fourier_transform sine_transform, inverse_sine_transform cosine_transform, inverse_cosine_transform hankel_transform, inverse_hankel_transform mellin_transform, laplace_transform """ return FourierTransform(f, x, k).doit(**hints) class InverseFourierTransform(FourierTypeTransform): """ Class representing unevaluated inverse Fourier transforms. For usage of this class, see the :class:`IntegralTransform` docstring. For how to compute inverse Fourier transforms, see the :func:`inverse_fourier_transform` docstring. """ _name = 'Inverse Fourier' def a(self): return 1 def b(self): return 2*S.Pi def inverse_fourier_transform(F, k, x, **hints): r""" Compute the unitary, ordinary-frequency inverse Fourier transform of `F`, defined as .. math:: f(x) = \int_{-\infty}^\infty F(k) e^{2\pi i x k} \mathrm{d} k. Explanation =========== If the transform cannot be computed in closed form, this function returns an unevaluated :class:`InverseFourierTransform` object. For other Fourier transform conventions, see the function :func:`sympy.integrals.transforms._fourier_transform`. For a description of possible hints, refer to the docstring of :func:`sympy.integrals.transforms.IntegralTransform.doit`. Note that for this transform, by default ``noconds=True``. Examples ======== >>> from sympy import inverse_fourier_transform, exp, sqrt, pi >>> from sympy.abc import x, k >>> inverse_fourier_transform(sqrt(pi)*exp(-(pi*k)**2), k, x) exp(-x**2) >>> inverse_fourier_transform(sqrt(pi)*exp(-(pi*k)**2), k, x, noconds=False) (exp(-x**2), True) See Also ======== fourier_transform sine_transform, inverse_sine_transform cosine_transform, inverse_cosine_transform hankel_transform, inverse_hankel_transform mellin_transform, laplace_transform """ return InverseFourierTransform(F, k, x).doit(**hints) ########################################################################## # Fourier Sine and Cosine Transform ########################################################################## @_noconds_(True) def _sine_cosine_transform(f, x, k, a, b, K, name, simplify=True): """ Compute a general sine or cosine-type transform F(k) = a int_0^oo b*sin(x*k) f(x) dx. F(k) = a int_0^oo b*cos(x*k) f(x) dx. For suitable choice of a and b, this reduces to the standard sine/cosine and inverse sine/cosine transforms. """ F = integrate(a*f*K(b*x*k), (x, S.Zero, S.Infinity)) if not F.has(Integral): return _simplify(F, simplify), S.true if not F.is_Piecewise: raise IntegralTransformError(name, f, 'could not compute integral') F, cond = F.args[0] if F.has(Integral): raise IntegralTransformError(name, f, 'integral in unexpected form') return _simplify(F, simplify), cond class SineCosineTypeTransform(IntegralTransform): """ Base class for sine and cosine transforms. Specify cls._kern. """ def a(self): raise NotImplementedError( "Class %s must implement a(self) but does not" % self.__class__) def b(self): raise NotImplementedError( "Class %s must implement b(self) but does not" % self.__class__) def _compute_transform(self, f, x, k, **hints): return _sine_cosine_transform(f, x, k, self.a(), self.b(), self.__class__._kern, self.__class__._name, **hints) def _as_integral(self, f, x, k): a = self.a() b = self.b() K = self.__class__._kern return Integral(a*f*K(b*x*k), (x, S.Zero, S.Infinity)) class SineTransform(SineCosineTypeTransform): """ Class representing unevaluated sine transforms. For usage of this class, see the :class:`IntegralTransform` docstring. For how to compute sine transforms, see the :func:`sine_transform` docstring. """ _name = 'Sine' _kern = sin def a(self): return sqrt(2)/sqrt(pi) def b(self): return S.One def sine_transform(f, x, k, **hints): r""" Compute the unitary, ordinary-frequency sine transform of `f`, defined as .. math:: F(k) = \sqrt{\frac{2}{\pi}} \int_{0}^\infty f(x) \sin(2\pi x k) \mathrm{d} x. Explanation =========== If the transform cannot be computed in closed form, this function returns an unevaluated :class:`SineTransform` object. For a description of possible hints, refer to the docstring of :func:`sympy.integrals.transforms.IntegralTransform.doit`. Note that for this transform, by default ``noconds=True``. Examples ======== >>> from sympy import sine_transform, exp >>> from sympy.abc import x, k, a >>> sine_transform(x*exp(-a*x**2), x, k) sqrt(2)*k*exp(-k**2/(4*a))/(4*a**(3/2)) >>> sine_transform(x**(-a), x, k) 2**(1/2 - a)*k**(a - 1)*gamma(1 - a/2)/gamma(a/2 + 1/2) See Also ======== fourier_transform, inverse_fourier_transform inverse_sine_transform cosine_transform, inverse_cosine_transform hankel_transform, inverse_hankel_transform mellin_transform, laplace_transform """ return SineTransform(f, x, k).doit(**hints) class InverseSineTransform(SineCosineTypeTransform): """ Class representing unevaluated inverse sine transforms. For usage of this class, see the :class:`IntegralTransform` docstring. For how to compute inverse sine transforms, see the :func:`inverse_sine_transform` docstring. """ _name = 'Inverse Sine' _kern = sin def a(self): return sqrt(2)/sqrt(pi) def b(self): return S.One def inverse_sine_transform(F, k, x, **hints): r""" Compute the unitary, ordinary-frequency inverse sine transform of `F`, defined as .. math:: f(x) = \sqrt{\frac{2}{\pi}} \int_{0}^\infty F(k) \sin(2\pi x k) \mathrm{d} k. Explanation =========== If the transform cannot be computed in closed form, this function returns an unevaluated :class:`InverseSineTransform` object. For a description of possible hints, refer to the docstring of :func:`sympy.integrals.transforms.IntegralTransform.doit`. Note that for this transform, by default ``noconds=True``. Examples ======== >>> from sympy import inverse_sine_transform, exp, sqrt, gamma >>> from sympy.abc import x, k, a >>> inverse_sine_transform(2**((1-2*a)/2)*k**(a - 1)* ... gamma(-a/2 + 1)/gamma((a+1)/2), k, x) x**(-a) >>> inverse_sine_transform(sqrt(2)*k*exp(-k**2/(4*a))/(4*sqrt(a)**3), k, x) x*exp(-a*x**2) See Also ======== fourier_transform, inverse_fourier_transform sine_transform cosine_transform, inverse_cosine_transform hankel_transform, inverse_hankel_transform mellin_transform, laplace_transform """ return InverseSineTransform(F, k, x).doit(**hints) class CosineTransform(SineCosineTypeTransform): """ Class representing unevaluated cosine transforms. For usage of this class, see the :class:`IntegralTransform` docstring. For how to compute cosine transforms, see the :func:`cosine_transform` docstring. """ _name = 'Cosine' _kern = cos def a(self): return sqrt(2)/sqrt(pi) def b(self): return S.One def cosine_transform(f, x, k, **hints): r""" Compute the unitary, ordinary-frequency cosine transform of `f`, defined as .. math:: F(k) = \sqrt{\frac{2}{\pi}} \int_{0}^\infty f(x) \cos(2\pi x k) \mathrm{d} x. Explanation =========== If the transform cannot be computed in closed form, this function returns an unevaluated :class:`CosineTransform` object. For a description of possible hints, refer to the docstring of :func:`sympy.integrals.transforms.IntegralTransform.doit`. Note that for this transform, by default ``noconds=True``. Examples ======== >>> from sympy import cosine_transform, exp, sqrt, cos >>> from sympy.abc import x, k, a >>> cosine_transform(exp(-a*x), x, k) sqrt(2)*a/(sqrt(pi)*(a**2 + k**2)) >>> cosine_transform(exp(-a*sqrt(x))*cos(a*sqrt(x)), x, k) a*exp(-a**2/(2*k))/(2*k**(3/2)) See Also ======== fourier_transform, inverse_fourier_transform, sine_transform, inverse_sine_transform inverse_cosine_transform hankel_transform, inverse_hankel_transform mellin_transform, laplace_transform """ return CosineTransform(f, x, k).doit(**hints) class InverseCosineTransform(SineCosineTypeTransform): """ Class representing unevaluated inverse cosine transforms. For usage of this class, see the :class:`IntegralTransform` docstring. For how to compute inverse cosine transforms, see the :func:`inverse_cosine_transform` docstring. """ _name = 'Inverse Cosine' _kern = cos def a(self): return sqrt(2)/sqrt(pi) def b(self): return S.One def inverse_cosine_transform(F, k, x, **hints): r""" Compute the unitary, ordinary-frequency inverse cosine transform of `F`, defined as .. math:: f(x) = \sqrt{\frac{2}{\pi}} \int_{0}^\infty F(k) \cos(2\pi x k) \mathrm{d} k. Explanation =========== If the transform cannot be computed in closed form, this function returns an unevaluated :class:`InverseCosineTransform` object. For a description of possible hints, refer to the docstring of :func:`sympy.integrals.transforms.IntegralTransform.doit`. Note that for this transform, by default ``noconds=True``. Examples ======== >>> from sympy import inverse_cosine_transform, sqrt, pi >>> from sympy.abc import x, k, a >>> inverse_cosine_transform(sqrt(2)*a/(sqrt(pi)*(a**2 + k**2)), k, x) exp(-a*x) >>> inverse_cosine_transform(1/sqrt(k), k, x) 1/sqrt(x) See Also ======== fourier_transform, inverse_fourier_transform, sine_transform, inverse_sine_transform cosine_transform hankel_transform, inverse_hankel_transform mellin_transform, laplace_transform """ return InverseCosineTransform(F, k, x).doit(**hints) ########################################################################## # Hankel Transform ########################################################################## @_noconds_(True) def _hankel_transform(f, r, k, nu, name, simplify=True): r""" Compute a general Hankel transform .. math:: F_\nu(k) = \int_{0}^\infty f(r) J_\nu(k r) r \mathrm{d} r. """ F = integrate(f*besselj(nu, k*r)*r, (r, S.Zero, S.Infinity)) if not F.has(Integral): return _simplify(F, simplify), S.true if not F.is_Piecewise: raise IntegralTransformError(name, f, 'could not compute integral') F, cond = F.args[0] if F.has(Integral): raise IntegralTransformError(name, f, 'integral in unexpected form') return _simplify(F, simplify), cond class HankelTypeTransform(IntegralTransform): """ Base class for Hankel transforms. """ def doit(self, **hints): return self._compute_transform(self.function, self.function_variable, self.transform_variable, self.args[3], **hints) def _compute_transform(self, f, r, k, nu, **hints): return _hankel_transform(f, r, k, nu, self._name, **hints) def _as_integral(self, f, r, k, nu): return Integral(f*besselj(nu, k*r)*r, (r, S.Zero, S.Infinity)) @property def as_integral(self): return self._as_integral(self.function, self.function_variable, self.transform_variable, self.args[3]) class HankelTransform(HankelTypeTransform): """ Class representing unevaluated Hankel transforms. For usage of this class, see the :class:`IntegralTransform` docstring. For how to compute Hankel transforms, see the :func:`hankel_transform` docstring. """ _name = 'Hankel' def hankel_transform(f, r, k, nu, **hints): r""" Compute the Hankel transform of `f`, defined as .. math:: F_\nu(k) = \int_{0}^\infty f(r) J_\nu(k r) r \mathrm{d} r. Explanation =========== If the transform cannot be computed in closed form, this function returns an unevaluated :class:`HankelTransform` object. For a description of possible hints, refer to the docstring of :func:`sympy.integrals.transforms.IntegralTransform.doit`. Note that for this transform, by default ``noconds=True``. Examples ======== >>> from sympy import hankel_transform, inverse_hankel_transform >>> from sympy import exp >>> from sympy.abc import r, k, m, nu, a >>> ht = hankel_transform(1/r**m, r, k, nu) >>> ht 2*k**(m - 2)*gamma(-m/2 + nu/2 + 1)/(2**m*gamma(m/2 + nu/2)) >>> inverse_hankel_transform(ht, k, r, nu) r**(-m) >>> ht = hankel_transform(exp(-a*r), r, k, 0) >>> ht a/(k**3*(a**2/k**2 + 1)**(3/2)) >>> inverse_hankel_transform(ht, k, r, 0) exp(-a*r) See Also ======== fourier_transform, inverse_fourier_transform sine_transform, inverse_sine_transform cosine_transform, inverse_cosine_transform inverse_hankel_transform mellin_transform, laplace_transform """ return HankelTransform(f, r, k, nu).doit(**hints) class InverseHankelTransform(HankelTypeTransform): """ Class representing unevaluated inverse Hankel transforms. For usage of this class, see the :class:`IntegralTransform` docstring. For how to compute inverse Hankel transforms, see the :func:`inverse_hankel_transform` docstring. """ _name = 'Inverse Hankel' def inverse_hankel_transform(F, k, r, nu, **hints): r""" Compute the inverse Hankel transform of `F` defined as .. math:: f(r) = \int_{0}^\infty F_\nu(k) J_\nu(k r) k \mathrm{d} k. Explanation =========== If the transform cannot be computed in closed form, this function returns an unevaluated :class:`InverseHankelTransform` object. For a description of possible hints, refer to the docstring of :func:`sympy.integrals.transforms.IntegralTransform.doit`. Note that for this transform, by default ``noconds=True``. Examples ======== >>> from sympy import hankel_transform, inverse_hankel_transform >>> from sympy import exp >>> from sympy.abc import r, k, m, nu, a >>> ht = hankel_transform(1/r**m, r, k, nu) >>> ht 2*k**(m - 2)*gamma(-m/2 + nu/2 + 1)/(2**m*gamma(m/2 + nu/2)) >>> inverse_hankel_transform(ht, k, r, nu) r**(-m) >>> ht = hankel_transform(exp(-a*r), r, k, 0) >>> ht a/(k**3*(a**2/k**2 + 1)**(3/2)) >>> inverse_hankel_transform(ht, k, r, 0) exp(-a*r) See Also ======== fourier_transform, inverse_fourier_transform sine_transform, inverse_sine_transform cosine_transform, inverse_cosine_transform hankel_transform mellin_transform, laplace_transform """ return InverseHankelTransform(F, k, r, nu).doit(**hints)