2191 lines
79 KiB
Python
2191 lines
79 KiB
Python
"""
|
|
Integrate functions by rewriting them as Meijer G-functions.
|
|
|
|
There are three user-visible functions that can be used by other parts of the
|
|
sympy library to solve various integration problems:
|
|
|
|
- meijerint_indefinite
|
|
- meijerint_definite
|
|
- meijerint_inversion
|
|
|
|
They can be used to compute, respectively, indefinite integrals, definite
|
|
integrals over intervals of the real line, and inverse laplace-type integrals
|
|
(from c-I*oo to c+I*oo). See the respective docstrings for details.
|
|
|
|
The main references for this are:
|
|
|
|
[L] Luke, Y. L. (1969), The Special Functions and Their Approximations,
|
|
Volume 1
|
|
|
|
[R] Kelly B. Roach. Meijer G Function Representations.
|
|
In: Proceedings of the 1997 International Symposium on Symbolic and
|
|
Algebraic Computation, pages 205-211, New York, 1997. ACM.
|
|
|
|
[P] A. P. Prudnikov, Yu. A. Brychkov and O. I. Marichev (1990).
|
|
Integrals and Series: More Special Functions, Vol. 3,.
|
|
Gordon and Breach Science Publisher
|
|
"""
|
|
|
|
from __future__ import annotations
|
|
import itertools
|
|
|
|
from sympy import SYMPY_DEBUG
|
|
from sympy.core import S, Expr
|
|
from sympy.core.add import Add
|
|
from sympy.core.basic import Basic
|
|
from sympy.core.cache import cacheit
|
|
from sympy.core.containers import Tuple
|
|
from sympy.core.exprtools import factor_terms
|
|
from sympy.core.function import (expand, expand_mul, expand_power_base,
|
|
expand_trig, Function)
|
|
from sympy.core.mul import Mul
|
|
from sympy.core.numbers import ilcm, Rational, pi
|
|
from sympy.core.relational import Eq, Ne, _canonical_coeff
|
|
from sympy.core.sorting import default_sort_key, ordered
|
|
from sympy.core.symbol import Dummy, symbols, Wild, Symbol
|
|
from sympy.core.sympify import sympify
|
|
from sympy.functions.combinatorial.factorials import factorial
|
|
from sympy.functions.elementary.complexes import (re, im, arg, Abs, sign,
|
|
unpolarify, polarify, polar_lift, principal_branch, unbranched_argument,
|
|
periodic_argument)
|
|
from sympy.functions.elementary.exponential import exp, exp_polar, log
|
|
from sympy.functions.elementary.integers import ceiling
|
|
from sympy.functions.elementary.hyperbolic import (cosh, sinh,
|
|
_rewrite_hyperbolics_as_exp, HyperbolicFunction)
|
|
from sympy.functions.elementary.miscellaneous import sqrt
|
|
from sympy.functions.elementary.piecewise import Piecewise, piecewise_fold
|
|
from sympy.functions.elementary.trigonometric import (cos, sin, sinc,
|
|
TrigonometricFunction)
|
|
from sympy.functions.special.bessel import besselj, bessely, besseli, besselk
|
|
from sympy.functions.special.delta_functions import DiracDelta, Heaviside
|
|
from sympy.functions.special.elliptic_integrals import elliptic_k, elliptic_e
|
|
from sympy.functions.special.error_functions import (erf, erfc, erfi, Ei,
|
|
expint, Si, Ci, Shi, Chi, fresnels, fresnelc)
|
|
from sympy.functions.special.gamma_functions import gamma
|
|
from sympy.functions.special.hyper import hyper, meijerg
|
|
from sympy.functions.special.singularity_functions import SingularityFunction
|
|
from .integrals import Integral
|
|
from sympy.logic.boolalg import And, Or, BooleanAtom, Not, BooleanFunction
|
|
from sympy.polys import cancel, factor
|
|
from sympy.utilities.iterables import multiset_partitions
|
|
from sympy.utilities.misc import debug as _debug
|
|
from sympy.utilities.misc import debugf as _debugf
|
|
|
|
# keep this at top for easy reference
|
|
z = Dummy('z')
|
|
|
|
|
|
def _has(res, *f):
|
|
# return True if res has f; in the case of Piecewise
|
|
# only return True if *all* pieces have f
|
|
res = piecewise_fold(res)
|
|
if getattr(res, 'is_Piecewise', False):
|
|
return all(_has(i, *f) for i in res.args)
|
|
return res.has(*f)
|
|
|
|
|
|
def _create_lookup_table(table):
|
|
""" Add formulae for the function -> meijerg lookup table. """
|
|
def wild(n):
|
|
return Wild(n, exclude=[z])
|
|
p, q, a, b, c = list(map(wild, 'pqabc'))
|
|
n = Wild('n', properties=[lambda x: x.is_Integer and x > 0])
|
|
t = p*z**q
|
|
|
|
def add(formula, an, ap, bm, bq, arg=t, fac=S.One, cond=True, hint=True):
|
|
table.setdefault(_mytype(formula, z), []).append((formula,
|
|
[(fac, meijerg(an, ap, bm, bq, arg))], cond, hint))
|
|
|
|
def addi(formula, inst, cond, hint=True):
|
|
table.setdefault(
|
|
_mytype(formula, z), []).append((formula, inst, cond, hint))
|
|
|
|
def constant(a):
|
|
return [(a, meijerg([1], [], [], [0], z)),
|
|
(a, meijerg([], [1], [0], [], z))]
|
|
table[()] = [(a, constant(a), True, True)]
|
|
|
|
# [P], Section 8.
|
|
class IsNonPositiveInteger(Function):
|
|
|
|
@classmethod
|
|
def eval(cls, arg):
|
|
arg = unpolarify(arg)
|
|
if arg.is_Integer is True:
|
|
return arg <= 0
|
|
|
|
# Section 8.4.2
|
|
# TODO this needs more polar_lift (c/f entry for exp)
|
|
add(Heaviside(t - b)*(t - b)**(a - 1), [a], [], [], [0], t/b,
|
|
gamma(a)*b**(a - 1), And(b > 0))
|
|
add(Heaviside(b - t)*(b - t)**(a - 1), [], [a], [0], [], t/b,
|
|
gamma(a)*b**(a - 1), And(b > 0))
|
|
add(Heaviside(z - (b/p)**(1/q))*(t - b)**(a - 1), [a], [], [], [0], t/b,
|
|
gamma(a)*b**(a - 1), And(b > 0))
|
|
add(Heaviside((b/p)**(1/q) - z)*(b - t)**(a - 1), [], [a], [0], [], t/b,
|
|
gamma(a)*b**(a - 1), And(b > 0))
|
|
add((b + t)**(-a), [1 - a], [], [0], [], t/b, b**(-a)/gamma(a),
|
|
hint=Not(IsNonPositiveInteger(a)))
|
|
add(Abs(b - t)**(-a), [1 - a], [(1 - a)/2], [0], [(1 - a)/2], t/b,
|
|
2*sin(pi*a/2)*gamma(1 - a)*Abs(b)**(-a), re(a) < 1)
|
|
add((t**a - b**a)/(t - b), [0, a], [], [0, a], [], t/b,
|
|
b**(a - 1)*sin(a*pi)/pi)
|
|
|
|
# 12
|
|
def A1(r, sign, nu):
|
|
return pi**Rational(-1, 2)*(-sign*nu/2)**(1 - 2*r)
|
|
|
|
def tmpadd(r, sgn):
|
|
# XXX the a**2 is bad for matching
|
|
add((sqrt(a**2 + t) + sgn*a)**b/(a**2 + t)**r,
|
|
[(1 + b)/2, 1 - 2*r + b/2], [],
|
|
[(b - sgn*b)/2], [(b + sgn*b)/2], t/a**2,
|
|
a**(b - 2*r)*A1(r, sgn, b))
|
|
tmpadd(0, 1)
|
|
tmpadd(0, -1)
|
|
tmpadd(S.Half, 1)
|
|
tmpadd(S.Half, -1)
|
|
|
|
# 13
|
|
def tmpadd(r, sgn):
|
|
add((sqrt(a + p*z**q) + sgn*sqrt(p)*z**(q/2))**b/(a + p*z**q)**r,
|
|
[1 - r + sgn*b/2], [1 - r - sgn*b/2], [0, S.Half], [],
|
|
p*z**q/a, a**(b/2 - r)*A1(r, sgn, b))
|
|
tmpadd(0, 1)
|
|
tmpadd(0, -1)
|
|
tmpadd(S.Half, 1)
|
|
tmpadd(S.Half, -1)
|
|
# (those after look obscure)
|
|
|
|
# Section 8.4.3
|
|
add(exp(polar_lift(-1)*t), [], [], [0], [])
|
|
|
|
# TODO can do sin^n, sinh^n by expansion ... where?
|
|
# 8.4.4 (hyperbolic functions)
|
|
add(sinh(t), [], [1], [S.Half], [1, 0], t**2/4, pi**Rational(3, 2))
|
|
add(cosh(t), [], [S.Half], [0], [S.Half, S.Half], t**2/4, pi**Rational(3, 2))
|
|
|
|
# Section 8.4.5
|
|
# TODO can do t + a. but can also do by expansion... (XXX not really)
|
|
add(sin(t), [], [], [S.Half], [0], t**2/4, sqrt(pi))
|
|
add(cos(t), [], [], [0], [S.Half], t**2/4, sqrt(pi))
|
|
|
|
# Section 8.4.6 (sinc function)
|
|
add(sinc(t), [], [], [0], [Rational(-1, 2)], t**2/4, sqrt(pi)/2)
|
|
|
|
# Section 8.5.5
|
|
def make_log1(subs):
|
|
N = subs[n]
|
|
return [(S.NegativeOne**N*factorial(N),
|
|
meijerg([], [1]*(N + 1), [0]*(N + 1), [], t))]
|
|
|
|
def make_log2(subs):
|
|
N = subs[n]
|
|
return [(factorial(N),
|
|
meijerg([1]*(N + 1), [], [], [0]*(N + 1), t))]
|
|
# TODO these only hold for positive p, and can be made more general
|
|
# but who uses log(x)*Heaviside(a-x) anyway ...
|
|
# TODO also it would be nice to derive them recursively ...
|
|
addi(log(t)**n*Heaviside(1 - t), make_log1, True)
|
|
addi(log(t)**n*Heaviside(t - 1), make_log2, True)
|
|
|
|
def make_log3(subs):
|
|
return make_log1(subs) + make_log2(subs)
|
|
addi(log(t)**n, make_log3, True)
|
|
addi(log(t + a),
|
|
constant(log(a)) + [(S.One, meijerg([1, 1], [], [1], [0], t/a))],
|
|
True)
|
|
addi(log(Abs(t - a)), constant(log(Abs(a))) +
|
|
[(pi, meijerg([1, 1], [S.Half], [1], [0, S.Half], t/a))],
|
|
True)
|
|
# TODO log(x)/(x+a) and log(x)/(x-1) can also be done. should they
|
|
# be derivable?
|
|
# TODO further formulae in this section seem obscure
|
|
|
|
# Sections 8.4.9-10
|
|
# TODO
|
|
|
|
# Section 8.4.11
|
|
addi(Ei(t),
|
|
constant(-S.ImaginaryUnit*pi) + [(S.NegativeOne, meijerg([], [1], [0, 0], [],
|
|
t*polar_lift(-1)))],
|
|
True)
|
|
|
|
# Section 8.4.12
|
|
add(Si(t), [1], [], [S.Half], [0, 0], t**2/4, sqrt(pi)/2)
|
|
add(Ci(t), [], [1], [0, 0], [S.Half], t**2/4, -sqrt(pi)/2)
|
|
|
|
# Section 8.4.13
|
|
add(Shi(t), [S.Half], [], [0], [Rational(-1, 2), Rational(-1, 2)], polar_lift(-1)*t**2/4,
|
|
t*sqrt(pi)/4)
|
|
add(Chi(t), [], [S.Half, 1], [0, 0], [S.Half, S.Half], t**2/4, -
|
|
pi**S('3/2')/2)
|
|
|
|
# generalized exponential integral
|
|
add(expint(a, t), [], [a], [a - 1, 0], [], t)
|
|
|
|
# Section 8.4.14
|
|
add(erf(t), [1], [], [S.Half], [0], t**2, 1/sqrt(pi))
|
|
# TODO exp(-x)*erf(I*x) does not work
|
|
add(erfc(t), [], [1], [0, S.Half], [], t**2, 1/sqrt(pi))
|
|
# This formula for erfi(z) yields a wrong(?) minus sign
|
|
#add(erfi(t), [1], [], [S.Half], [0], -t**2, I/sqrt(pi))
|
|
add(erfi(t), [S.Half], [], [0], [Rational(-1, 2)], -t**2, t/sqrt(pi))
|
|
|
|
# Fresnel Integrals
|
|
add(fresnels(t), [1], [], [Rational(3, 4)], [0, Rational(1, 4)], pi**2*t**4/16, S.Half)
|
|
add(fresnelc(t), [1], [], [Rational(1, 4)], [0, Rational(3, 4)], pi**2*t**4/16, S.Half)
|
|
|
|
##### bessel-type functions #####
|
|
# Section 8.4.19
|
|
add(besselj(a, t), [], [], [a/2], [-a/2], t**2/4)
|
|
|
|
# all of the following are derivable
|
|
#add(sin(t)*besselj(a, t), [Rational(1, 4), Rational(3, 4)], [], [(1+a)/2],
|
|
# [-a/2, a/2, (1-a)/2], t**2, 1/sqrt(2))
|
|
#add(cos(t)*besselj(a, t), [Rational(1, 4), Rational(3, 4)], [], [a/2],
|
|
# [-a/2, (1+a)/2, (1-a)/2], t**2, 1/sqrt(2))
|
|
#add(besselj(a, t)**2, [S.Half], [], [a], [-a, 0], t**2, 1/sqrt(pi))
|
|
#add(besselj(a, t)*besselj(b, t), [0, S.Half], [], [(a + b)/2],
|
|
# [-(a+b)/2, (a - b)/2, (b - a)/2], t**2, 1/sqrt(pi))
|
|
|
|
# Section 8.4.20
|
|
add(bessely(a, t), [], [-(a + 1)/2], [a/2, -a/2], [-(a + 1)/2], t**2/4)
|
|
|
|
# TODO all of the following should be derivable
|
|
#add(sin(t)*bessely(a, t), [Rational(1, 4), Rational(3, 4)], [(1 - a - 1)/2],
|
|
# [(1 + a)/2, (1 - a)/2], [(1 - a - 1)/2, (1 - 1 - a)/2, (1 - 1 + a)/2],
|
|
# t**2, 1/sqrt(2))
|
|
#add(cos(t)*bessely(a, t), [Rational(1, 4), Rational(3, 4)], [(0 - a - 1)/2],
|
|
# [(0 + a)/2, (0 - a)/2], [(0 - a - 1)/2, (1 - 0 - a)/2, (1 - 0 + a)/2],
|
|
# t**2, 1/sqrt(2))
|
|
#add(besselj(a, t)*bessely(b, t), [0, S.Half], [(a - b - 1)/2],
|
|
# [(a + b)/2, (a - b)/2], [(a - b - 1)/2, -(a + b)/2, (b - a)/2],
|
|
# t**2, 1/sqrt(pi))
|
|
#addi(bessely(a, t)**2,
|
|
# [(2/sqrt(pi), meijerg([], [S.Half, S.Half - a], [0, a, -a],
|
|
# [S.Half - a], t**2)),
|
|
# (1/sqrt(pi), meijerg([S.Half], [], [a], [-a, 0], t**2))],
|
|
# True)
|
|
#addi(bessely(a, t)*bessely(b, t),
|
|
# [(2/sqrt(pi), meijerg([], [0, S.Half, (1 - a - b)/2],
|
|
# [(a + b)/2, (a - b)/2, (b - a)/2, -(a + b)/2],
|
|
# [(1 - a - b)/2], t**2)),
|
|
# (1/sqrt(pi), meijerg([0, S.Half], [], [(a + b)/2],
|
|
# [-(a + b)/2, (a - b)/2, (b - a)/2], t**2))],
|
|
# True)
|
|
|
|
# Section 8.4.21 ?
|
|
# Section 8.4.22
|
|
add(besseli(a, t), [], [(1 + a)/2], [a/2], [-a/2, (1 + a)/2], t**2/4, pi)
|
|
# TODO many more formulas. should all be derivable
|
|
|
|
# Section 8.4.23
|
|
add(besselk(a, t), [], [], [a/2, -a/2], [], t**2/4, S.Half)
|
|
# TODO many more formulas. should all be derivable
|
|
|
|
# Complete elliptic integrals K(z) and E(z)
|
|
add(elliptic_k(t), [S.Half, S.Half], [], [0], [0], -t, S.Half)
|
|
add(elliptic_e(t), [S.Half, 3*S.Half], [], [0], [0], -t, Rational(-1, 2)/2)
|
|
|
|
|
|
####################################################################
|
|
# First some helper functions.
|
|
####################################################################
|
|
|
|
from sympy.utilities.timeutils import timethis
|
|
timeit = timethis('meijerg')
|
|
|
|
|
|
def _mytype(f: Basic, x: Symbol) -> tuple[type[Basic], ...]:
|
|
""" Create a hashable entity describing the type of f. """
|
|
def key(x: type[Basic]) -> tuple[int, int, str]:
|
|
return x.class_key()
|
|
|
|
if x not in f.free_symbols:
|
|
return ()
|
|
elif f.is_Function:
|
|
return type(f),
|
|
return tuple(sorted((t for a in f.args for t in _mytype(a, x)), key=key))
|
|
|
|
|
|
class _CoeffExpValueError(ValueError):
|
|
"""
|
|
Exception raised by _get_coeff_exp, for internal use only.
|
|
"""
|
|
pass
|
|
|
|
|
|
def _get_coeff_exp(expr, x):
|
|
"""
|
|
When expr is known to be of the form c*x**b, with c and/or b possibly 1,
|
|
return c, b.
|
|
|
|
Examples
|
|
========
|
|
|
|
>>> from sympy.abc import x, a, b
|
|
>>> from sympy.integrals.meijerint import _get_coeff_exp
|
|
>>> _get_coeff_exp(a*x**b, x)
|
|
(a, b)
|
|
>>> _get_coeff_exp(x, x)
|
|
(1, 1)
|
|
>>> _get_coeff_exp(2*x, x)
|
|
(2, 1)
|
|
>>> _get_coeff_exp(x**3, x)
|
|
(1, 3)
|
|
"""
|
|
from sympy.simplify import powsimp
|
|
(c, m) = expand_power_base(powsimp(expr)).as_coeff_mul(x)
|
|
if not m:
|
|
return c, S.Zero
|
|
[m] = m
|
|
if m.is_Pow:
|
|
if m.base != x:
|
|
raise _CoeffExpValueError('expr not of form a*x**b')
|
|
return c, m.exp
|
|
elif m == x:
|
|
return c, S.One
|
|
else:
|
|
raise _CoeffExpValueError('expr not of form a*x**b: %s' % expr)
|
|
|
|
|
|
def _exponents(expr, x):
|
|
"""
|
|
Find the exponents of ``x`` (not including zero) in ``expr``.
|
|
|
|
Examples
|
|
========
|
|
|
|
>>> from sympy.integrals.meijerint import _exponents
|
|
>>> from sympy.abc import x, y
|
|
>>> from sympy import sin
|
|
>>> _exponents(x, x)
|
|
{1}
|
|
>>> _exponents(x**2, x)
|
|
{2}
|
|
>>> _exponents(x**2 + x, x)
|
|
{1, 2}
|
|
>>> _exponents(x**3*sin(x + x**y) + 1/x, x)
|
|
{-1, 1, 3, y}
|
|
"""
|
|
def _exponents_(expr, x, res):
|
|
if expr == x:
|
|
res.update([1])
|
|
return
|
|
if expr.is_Pow and expr.base == x:
|
|
res.update([expr.exp])
|
|
return
|
|
for argument in expr.args:
|
|
_exponents_(argument, x, res)
|
|
res = set()
|
|
_exponents_(expr, x, res)
|
|
return res
|
|
|
|
|
|
def _functions(expr, x):
|
|
""" Find the types of functions in expr, to estimate the complexity. """
|
|
return {e.func for e in expr.atoms(Function) if x in e.free_symbols}
|
|
|
|
|
|
def _find_splitting_points(expr, x):
|
|
"""
|
|
Find numbers a such that a linear substitution x -> x + a would
|
|
(hopefully) simplify expr.
|
|
|
|
Examples
|
|
========
|
|
|
|
>>> from sympy.integrals.meijerint import _find_splitting_points as fsp
|
|
>>> from sympy import sin
|
|
>>> from sympy.abc import x
|
|
>>> fsp(x, x)
|
|
{0}
|
|
>>> fsp((x-1)**3, x)
|
|
{1}
|
|
>>> fsp(sin(x+3)*x, x)
|
|
{-3, 0}
|
|
"""
|
|
p, q = [Wild(n, exclude=[x]) for n in 'pq']
|
|
|
|
def compute_innermost(expr, res):
|
|
if not isinstance(expr, Expr):
|
|
return
|
|
m = expr.match(p*x + q)
|
|
if m and m[p] != 0:
|
|
res.add(-m[q]/m[p])
|
|
return
|
|
if expr.is_Atom:
|
|
return
|
|
for argument in expr.args:
|
|
compute_innermost(argument, res)
|
|
innermost = set()
|
|
compute_innermost(expr, innermost)
|
|
return innermost
|
|
|
|
|
|
def _split_mul(f, x):
|
|
"""
|
|
Split expression ``f`` into fac, po, g, where fac is a constant factor,
|
|
po = x**s for some s independent of s, and g is "the rest".
|
|
|
|
Examples
|
|
========
|
|
|
|
>>> from sympy.integrals.meijerint import _split_mul
|
|
>>> from sympy import sin
|
|
>>> from sympy.abc import s, x
|
|
>>> _split_mul((3*x)**s*sin(x**2)*x, x)
|
|
(3**s, x*x**s, sin(x**2))
|
|
"""
|
|
fac = S.One
|
|
po = S.One
|
|
g = S.One
|
|
f = expand_power_base(f)
|
|
|
|
args = Mul.make_args(f)
|
|
for a in args:
|
|
if a == x:
|
|
po *= x
|
|
elif x not in a.free_symbols:
|
|
fac *= a
|
|
else:
|
|
if a.is_Pow and x not in a.exp.free_symbols:
|
|
c, t = a.base.as_coeff_mul(x)
|
|
if t != (x,):
|
|
c, t = expand_mul(a.base).as_coeff_mul(x)
|
|
if t == (x,):
|
|
po *= x**a.exp
|
|
fac *= unpolarify(polarify(c**a.exp, subs=False))
|
|
continue
|
|
g *= a
|
|
|
|
return fac, po, g
|
|
|
|
|
|
def _mul_args(f):
|
|
"""
|
|
Return a list ``L`` such that ``Mul(*L) == f``.
|
|
|
|
If ``f`` is not a ``Mul`` or ``Pow``, ``L=[f]``.
|
|
If ``f=g**n`` for an integer ``n``, ``L=[g]*n``.
|
|
If ``f`` is a ``Mul``, ``L`` comes from applying ``_mul_args`` to all factors of ``f``.
|
|
"""
|
|
args = Mul.make_args(f)
|
|
gs = []
|
|
for g in args:
|
|
if g.is_Pow and g.exp.is_Integer:
|
|
n = g.exp
|
|
base = g.base
|
|
if n < 0:
|
|
n = -n
|
|
base = 1/base
|
|
gs += [base]*n
|
|
else:
|
|
gs.append(g)
|
|
return gs
|
|
|
|
|
|
def _mul_as_two_parts(f):
|
|
"""
|
|
Find all the ways to split ``f`` into a product of two terms.
|
|
Return None on failure.
|
|
|
|
Explanation
|
|
===========
|
|
|
|
Although the order is canonical from multiset_partitions, this is
|
|
not necessarily the best order to process the terms. For example,
|
|
if the case of len(gs) == 2 is removed and multiset is allowed to
|
|
sort the terms, some tests fail.
|
|
|
|
Examples
|
|
========
|
|
|
|
>>> from sympy.integrals.meijerint import _mul_as_two_parts
|
|
>>> from sympy import sin, exp, ordered
|
|
>>> from sympy.abc import x
|
|
>>> list(ordered(_mul_as_two_parts(x*sin(x)*exp(x))))
|
|
[(x, exp(x)*sin(x)), (x*exp(x), sin(x)), (x*sin(x), exp(x))]
|
|
"""
|
|
|
|
gs = _mul_args(f)
|
|
if len(gs) < 2:
|
|
return None
|
|
if len(gs) == 2:
|
|
return [tuple(gs)]
|
|
return [(Mul(*x), Mul(*y)) for (x, y) in multiset_partitions(gs, 2)]
|
|
|
|
|
|
def _inflate_g(g, n):
|
|
""" Return C, h such that h is a G function of argument z**n and
|
|
g = C*h. """
|
|
# TODO should this be a method of meijerg?
|
|
# See: [L, page 150, equation (5)]
|
|
def inflate(params, n):
|
|
""" (a1, .., ak) -> (a1/n, (a1+1)/n, ..., (ak + n-1)/n) """
|
|
return [(a + i)/n for a, i in itertools.product(params, range(n))]
|
|
v = S(len(g.ap) - len(g.bq))
|
|
C = n**(1 + g.nu + v/2)
|
|
C /= (2*pi)**((n - 1)*g.delta)
|
|
return C, meijerg(inflate(g.an, n), inflate(g.aother, n),
|
|
inflate(g.bm, n), inflate(g.bother, n),
|
|
g.argument**n * n**(n*v))
|
|
|
|
|
|
def _flip_g(g):
|
|
""" Turn the G function into one of inverse argument
|
|
(i.e. G(1/x) -> G'(x)) """
|
|
# See [L], section 5.2
|
|
def tr(l):
|
|
return [1 - a for a in l]
|
|
return meijerg(tr(g.bm), tr(g.bother), tr(g.an), tr(g.aother), 1/g.argument)
|
|
|
|
|
|
def _inflate_fox_h(g, a):
|
|
r"""
|
|
Let d denote the integrand in the definition of the G function ``g``.
|
|
Consider the function H which is defined in the same way, but with
|
|
integrand d/Gamma(a*s) (contour conventions as usual).
|
|
|
|
If ``a`` is rational, the function H can be written as C*G, for a constant C
|
|
and a G-function G.
|
|
|
|
This function returns C, G.
|
|
"""
|
|
if a < 0:
|
|
return _inflate_fox_h(_flip_g(g), -a)
|
|
p = S(a.p)
|
|
q = S(a.q)
|
|
# We use the substitution s->qs, i.e. inflate g by q. We are left with an
|
|
# extra factor of Gamma(p*s), for which we use Gauss' multiplication
|
|
# theorem.
|
|
D, g = _inflate_g(g, q)
|
|
z = g.argument
|
|
D /= (2*pi)**((1 - p)/2)*p**Rational(-1, 2)
|
|
z /= p**p
|
|
bs = [(n + 1)/p for n in range(p)]
|
|
return D, meijerg(g.an, g.aother, g.bm, list(g.bother) + bs, z)
|
|
|
|
_dummies: dict[tuple[str, str], Dummy] = {}
|
|
|
|
|
|
def _dummy(name, token, expr, **kwargs):
|
|
"""
|
|
Return a dummy. This will return the same dummy if the same token+name is
|
|
requested more than once, and it is not already in expr.
|
|
This is for being cache-friendly.
|
|
"""
|
|
d = _dummy_(name, token, **kwargs)
|
|
if d in expr.free_symbols:
|
|
return Dummy(name, **kwargs)
|
|
return d
|
|
|
|
|
|
def _dummy_(name, token, **kwargs):
|
|
"""
|
|
Return a dummy associated to name and token. Same effect as declaring
|
|
it globally.
|
|
"""
|
|
global _dummies
|
|
if not (name, token) in _dummies:
|
|
_dummies[(name, token)] = Dummy(name, **kwargs)
|
|
return _dummies[(name, token)]
|
|
|
|
|
|
def _is_analytic(f, x):
|
|
""" Check if f(x), when expressed using G functions on the positive reals,
|
|
will in fact agree with the G functions almost everywhere """
|
|
return not any(x in expr.free_symbols for expr in f.atoms(Heaviside, Abs))
|
|
|
|
|
|
def _condsimp(cond, first=True):
|
|
"""
|
|
Do naive simplifications on ``cond``.
|
|
|
|
Explanation
|
|
===========
|
|
|
|
Note that this routine is completely ad-hoc, simplification rules being
|
|
added as need arises rather than following any logical pattern.
|
|
|
|
Examples
|
|
========
|
|
|
|
>>> from sympy.integrals.meijerint import _condsimp as simp
|
|
>>> from sympy import Or, Eq
|
|
>>> from sympy.abc import x, y
|
|
>>> simp(Or(x < y, Eq(x, y)))
|
|
x <= y
|
|
"""
|
|
if first:
|
|
cond = cond.replace(lambda _: _.is_Relational, _canonical_coeff)
|
|
first = False
|
|
if not isinstance(cond, BooleanFunction):
|
|
return cond
|
|
p, q, r = symbols('p q r', cls=Wild)
|
|
# transforms tests use 0, 4, 5 and 11-14
|
|
# meijer tests use 0, 2, 11, 14
|
|
# joint_rv uses 6, 7
|
|
rules = [
|
|
(Or(p < q, Eq(p, q)), p <= q), # 0
|
|
# The next two obviously are instances of a general pattern, but it is
|
|
# easier to spell out the few cases we care about.
|
|
(And(Abs(arg(p)) <= pi, Abs(arg(p) - 2*pi) <= pi),
|
|
Eq(arg(p) - pi, 0)), # 1
|
|
(And(Abs(2*arg(p) + pi) <= pi, Abs(2*arg(p) - pi) <= pi),
|
|
Eq(arg(p), 0)), # 2
|
|
(And(Abs(2*arg(p) + pi) < pi, Abs(2*arg(p) - pi) <= pi),
|
|
S.false), # 3
|
|
(And(Abs(arg(p) - pi/2) <= pi/2, Abs(arg(p) + pi/2) <= pi/2),
|
|
Eq(arg(p), 0)), # 4
|
|
(And(Abs(arg(p) - pi/2) <= pi/2, Abs(arg(p) + pi/2) < pi/2),
|
|
S.false), # 5
|
|
(And(Abs(arg(p**2/2 + 1)) < pi, Ne(Abs(arg(p**2/2 + 1)), pi)),
|
|
S.true), # 6
|
|
(Or(Abs(arg(p**2/2 + 1)) < pi, Ne(1/(p**2/2 + 1), 0)),
|
|
S.true), # 7
|
|
(And(Abs(unbranched_argument(p)) <= pi,
|
|
Abs(unbranched_argument(exp_polar(-2*pi*S.ImaginaryUnit)*p)) <= pi),
|
|
Eq(unbranched_argument(exp_polar(-S.ImaginaryUnit*pi)*p), 0)), # 8
|
|
(And(Abs(unbranched_argument(p)) <= pi/2,
|
|
Abs(unbranched_argument(exp_polar(-pi*S.ImaginaryUnit)*p)) <= pi/2),
|
|
Eq(unbranched_argument(exp_polar(-S.ImaginaryUnit*pi/2)*p), 0)), # 9
|
|
(Or(p <= q, And(p < q, r)), p <= q), # 10
|
|
(Ne(p**2, 1) & (p**2 > 1), p**2 > 1), # 11
|
|
(Ne(1/p, 1) & (cos(Abs(arg(p)))*Abs(p) > 1), Abs(p) > 1), # 12
|
|
(Ne(p, 2) & (cos(Abs(arg(p)))*Abs(p) > 2), Abs(p) > 2), # 13
|
|
((Abs(arg(p)) < pi/2) & (cos(Abs(arg(p)))*sqrt(Abs(p**2)) > 1), p**2 > 1), # 14
|
|
]
|
|
cond = cond.func(*[_condsimp(_, first) for _ in cond.args])
|
|
change = True
|
|
while change:
|
|
change = False
|
|
for irule, (fro, to) in enumerate(rules):
|
|
if fro.func != cond.func:
|
|
continue
|
|
for n, arg1 in enumerate(cond.args):
|
|
if r in fro.args[0].free_symbols:
|
|
m = arg1.match(fro.args[1])
|
|
num = 1
|
|
else:
|
|
num = 0
|
|
m = arg1.match(fro.args[0])
|
|
if not m:
|
|
continue
|
|
otherargs = [x.subs(m) for x in fro.args[:num] + fro.args[num + 1:]]
|
|
otherlist = [n]
|
|
for arg2 in otherargs:
|
|
for k, arg3 in enumerate(cond.args):
|
|
if k in otherlist:
|
|
continue
|
|
if arg2 == arg3:
|
|
otherlist += [k]
|
|
break
|
|
if isinstance(arg3, And) and arg2.args[1] == r and \
|
|
isinstance(arg2, And) and arg2.args[0] in arg3.args:
|
|
otherlist += [k]
|
|
break
|
|
if isinstance(arg3, And) and arg2.args[0] == r and \
|
|
isinstance(arg2, And) and arg2.args[1] in arg3.args:
|
|
otherlist += [k]
|
|
break
|
|
if len(otherlist) != len(otherargs) + 1:
|
|
continue
|
|
newargs = [arg_ for (k, arg_) in enumerate(cond.args)
|
|
if k not in otherlist] + [to.subs(m)]
|
|
if SYMPY_DEBUG:
|
|
if irule not in (0, 2, 4, 5, 6, 7, 11, 12, 13, 14):
|
|
print('used new rule:', irule)
|
|
cond = cond.func(*newargs)
|
|
change = True
|
|
break
|
|
|
|
# final tweak
|
|
def rel_touchup(rel):
|
|
if rel.rel_op != '==' or rel.rhs != 0:
|
|
return rel
|
|
|
|
# handle Eq(*, 0)
|
|
LHS = rel.lhs
|
|
m = LHS.match(arg(p)**q)
|
|
if not m:
|
|
m = LHS.match(unbranched_argument(polar_lift(p)**q))
|
|
if not m:
|
|
if isinstance(LHS, periodic_argument) and not LHS.args[0].is_polar \
|
|
and LHS.args[1] is S.Infinity:
|
|
return (LHS.args[0] > 0)
|
|
return rel
|
|
return (m[p] > 0)
|
|
cond = cond.replace(lambda _: _.is_Relational, rel_touchup)
|
|
if SYMPY_DEBUG:
|
|
print('_condsimp: ', cond)
|
|
return cond
|
|
|
|
def _eval_cond(cond):
|
|
""" Re-evaluate the conditions. """
|
|
if isinstance(cond, bool):
|
|
return cond
|
|
return _condsimp(cond.doit())
|
|
|
|
####################################################################
|
|
# Now the "backbone" functions to do actual integration.
|
|
####################################################################
|
|
|
|
|
|
def _my_principal_branch(expr, period, full_pb=False):
|
|
""" Bring expr nearer to its principal branch by removing superfluous
|
|
factors.
|
|
This function does *not* guarantee to yield the principal branch,
|
|
to avoid introducing opaque principal_branch() objects,
|
|
unless full_pb=True. """
|
|
res = principal_branch(expr, period)
|
|
if not full_pb:
|
|
res = res.replace(principal_branch, lambda x, y: x)
|
|
return res
|
|
|
|
|
|
def _rewrite_saxena_1(fac, po, g, x):
|
|
"""
|
|
Rewrite the integral fac*po*g dx, from zero to infinity, as
|
|
integral fac*G, where G has argument a*x. Note po=x**s.
|
|
Return fac, G.
|
|
"""
|
|
_, s = _get_coeff_exp(po, x)
|
|
a, b = _get_coeff_exp(g.argument, x)
|
|
period = g.get_period()
|
|
a = _my_principal_branch(a, period)
|
|
|
|
# We substitute t = x**b.
|
|
C = fac/(Abs(b)*a**((s + 1)/b - 1))
|
|
# Absorb a factor of (at)**((1 + s)/b - 1).
|
|
|
|
def tr(l):
|
|
return [a + (1 + s)/b - 1 for a in l]
|
|
return C, meijerg(tr(g.an), tr(g.aother), tr(g.bm), tr(g.bother),
|
|
a*x)
|
|
|
|
|
|
def _check_antecedents_1(g, x, helper=False):
|
|
r"""
|
|
Return a condition under which the mellin transform of g exists.
|
|
Any power of x has already been absorbed into the G function,
|
|
so this is just $\int_0^\infty g\, dx$.
|
|
|
|
See [L, section 5.6.1]. (Note that s=1.)
|
|
|
|
If ``helper`` is True, only check if the MT exists at infinity, i.e. if
|
|
$\int_1^\infty g\, dx$ exists.
|
|
"""
|
|
# NOTE if you update these conditions, please update the documentation as well
|
|
delta = g.delta
|
|
eta, _ = _get_coeff_exp(g.argument, x)
|
|
m, n, p, q = S([len(g.bm), len(g.an), len(g.ap), len(g.bq)])
|
|
|
|
if p > q:
|
|
def tr(l):
|
|
return [1 - x for x in l]
|
|
return _check_antecedents_1(meijerg(tr(g.bm), tr(g.bother),
|
|
tr(g.an), tr(g.aother), x/eta),
|
|
x)
|
|
|
|
tmp = [-re(b) < 1 for b in g.bm] + [1 < 1 - re(a) for a in g.an]
|
|
cond_3 = And(*tmp)
|
|
|
|
tmp += [-re(b) < 1 for b in g.bother]
|
|
tmp += [1 < 1 - re(a) for a in g.aother]
|
|
cond_3_star = And(*tmp)
|
|
|
|
cond_4 = (-re(g.nu) + (q + 1 - p)/2 > q - p)
|
|
|
|
def debug(*msg):
|
|
_debug(*msg)
|
|
|
|
def debugf(string, arg):
|
|
_debugf(string, arg)
|
|
|
|
debug('Checking antecedents for 1 function:')
|
|
debugf(' delta=%s, eta=%s, m=%s, n=%s, p=%s, q=%s',
|
|
(delta, eta, m, n, p, q))
|
|
debugf(' ap = %s, %s', (list(g.an), list(g.aother)))
|
|
debugf(' bq = %s, %s', (list(g.bm), list(g.bother)))
|
|
debugf(' cond_3=%s, cond_3*=%s, cond_4=%s', (cond_3, cond_3_star, cond_4))
|
|
|
|
conds = []
|
|
|
|
# case 1
|
|
case1 = []
|
|
tmp1 = [1 <= n, p < q, 1 <= m]
|
|
tmp2 = [1 <= p, 1 <= m, Eq(q, p + 1), Not(And(Eq(n, 0), Eq(m, p + 1)))]
|
|
tmp3 = [1 <= p, Eq(q, p)]
|
|
for k in range(ceiling(delta/2) + 1):
|
|
tmp3 += [Ne(Abs(unbranched_argument(eta)), (delta - 2*k)*pi)]
|
|
tmp = [delta > 0, Abs(unbranched_argument(eta)) < delta*pi]
|
|
extra = [Ne(eta, 0), cond_3]
|
|
if helper:
|
|
extra = []
|
|
for t in [tmp1, tmp2, tmp3]:
|
|
case1 += [And(*(t + tmp + extra))]
|
|
conds += case1
|
|
debug(' case 1:', case1)
|
|
|
|
# case 2
|
|
extra = [cond_3]
|
|
if helper:
|
|
extra = []
|
|
case2 = [And(Eq(n, 0), p + 1 <= m, m <= q,
|
|
Abs(unbranched_argument(eta)) < delta*pi, *extra)]
|
|
conds += case2
|
|
debug(' case 2:', case2)
|
|
|
|
# case 3
|
|
extra = [cond_3, cond_4]
|
|
if helper:
|
|
extra = []
|
|
case3 = [And(p < q, 1 <= m, delta > 0, Eq(Abs(unbranched_argument(eta)), delta*pi),
|
|
*extra)]
|
|
case3 += [And(p <= q - 2, Eq(delta, 0), Eq(Abs(unbranched_argument(eta)), 0), *extra)]
|
|
conds += case3
|
|
debug(' case 3:', case3)
|
|
|
|
# TODO altered cases 4-7
|
|
|
|
# extra case from wofram functions site:
|
|
# (reproduced verbatim from Prudnikov, section 2.24.2)
|
|
# https://functions.wolfram.com/HypergeometricFunctions/MeijerG/21/02/01/
|
|
case_extra = []
|
|
case_extra += [Eq(p, q), Eq(delta, 0), Eq(unbranched_argument(eta), 0), Ne(eta, 0)]
|
|
if not helper:
|
|
case_extra += [cond_3]
|
|
s = []
|
|
for a, b in zip(g.ap, g.bq):
|
|
s += [b - a]
|
|
case_extra += [re(Add(*s)) < 0]
|
|
case_extra = And(*case_extra)
|
|
conds += [case_extra]
|
|
debug(' extra case:', [case_extra])
|
|
|
|
case_extra_2 = [And(delta > 0, Abs(unbranched_argument(eta)) < delta*pi)]
|
|
if not helper:
|
|
case_extra_2 += [cond_3]
|
|
case_extra_2 = And(*case_extra_2)
|
|
conds += [case_extra_2]
|
|
debug(' second extra case:', [case_extra_2])
|
|
|
|
# TODO This leaves only one case from the three listed by Prudnikov.
|
|
# Investigate if these indeed cover everything; if so, remove the rest.
|
|
|
|
return Or(*conds)
|
|
|
|
|
|
def _int0oo_1(g, x):
|
|
r"""
|
|
Evaluate $\int_0^\infty g\, dx$ using G functions,
|
|
assuming the necessary conditions are fulfilled.
|
|
|
|
Examples
|
|
========
|
|
|
|
>>> from sympy.abc import a, b, c, d, x, y
|
|
>>> from sympy import meijerg
|
|
>>> from sympy.integrals.meijerint import _int0oo_1
|
|
>>> _int0oo_1(meijerg([a], [b], [c], [d], x*y), x)
|
|
gamma(-a)*gamma(c + 1)/(y*gamma(-d)*gamma(b + 1))
|
|
"""
|
|
from sympy.simplify import gammasimp
|
|
# See [L, section 5.6.1]. Note that s=1.
|
|
eta, _ = _get_coeff_exp(g.argument, x)
|
|
res = 1/eta
|
|
# XXX TODO we should reduce order first
|
|
for b in g.bm:
|
|
res *= gamma(b + 1)
|
|
for a in g.an:
|
|
res *= gamma(1 - a - 1)
|
|
for b in g.bother:
|
|
res /= gamma(1 - b - 1)
|
|
for a in g.aother:
|
|
res /= gamma(a + 1)
|
|
return gammasimp(unpolarify(res))
|
|
|
|
|
|
def _rewrite_saxena(fac, po, g1, g2, x, full_pb=False):
|
|
"""
|
|
Rewrite the integral ``fac*po*g1*g2`` from 0 to oo in terms of G
|
|
functions with argument ``c*x``.
|
|
|
|
Explanation
|
|
===========
|
|
|
|
Return C, f1, f2 such that integral C f1 f2 from 0 to infinity equals
|
|
integral fac ``po``, ``g1``, ``g2`` from 0 to infinity.
|
|
|
|
Examples
|
|
========
|
|
|
|
>>> from sympy.integrals.meijerint import _rewrite_saxena
|
|
>>> from sympy.abc import s, t, m
|
|
>>> from sympy import meijerg
|
|
>>> g1 = meijerg([], [], [0], [], s*t)
|
|
>>> g2 = meijerg([], [], [m/2], [-m/2], t**2/4)
|
|
>>> r = _rewrite_saxena(1, t**0, g1, g2, t)
|
|
>>> r[0]
|
|
s/(4*sqrt(pi))
|
|
>>> r[1]
|
|
meijerg(((), ()), ((-1/2, 0), ()), s**2*t/4)
|
|
>>> r[2]
|
|
meijerg(((), ()), ((m/2,), (-m/2,)), t/4)
|
|
"""
|
|
def pb(g):
|
|
a, b = _get_coeff_exp(g.argument, x)
|
|
per = g.get_period()
|
|
return meijerg(g.an, g.aother, g.bm, g.bother,
|
|
_my_principal_branch(a, per, full_pb)*x**b)
|
|
|
|
_, s = _get_coeff_exp(po, x)
|
|
_, b1 = _get_coeff_exp(g1.argument, x)
|
|
_, b2 = _get_coeff_exp(g2.argument, x)
|
|
if (b1 < 0) == True:
|
|
b1 = -b1
|
|
g1 = _flip_g(g1)
|
|
if (b2 < 0) == True:
|
|
b2 = -b2
|
|
g2 = _flip_g(g2)
|
|
if not b1.is_Rational or not b2.is_Rational:
|
|
return
|
|
m1, n1 = b1.p, b1.q
|
|
m2, n2 = b2.p, b2.q
|
|
tau = ilcm(m1*n2, m2*n1)
|
|
r1 = tau//(m1*n2)
|
|
r2 = tau//(m2*n1)
|
|
|
|
C1, g1 = _inflate_g(g1, r1)
|
|
C2, g2 = _inflate_g(g2, r2)
|
|
g1 = pb(g1)
|
|
g2 = pb(g2)
|
|
|
|
fac *= C1*C2
|
|
a1, b = _get_coeff_exp(g1.argument, x)
|
|
a2, _ = _get_coeff_exp(g2.argument, x)
|
|
|
|
# arbitrarily tack on the x**s part to g1
|
|
# TODO should we try both?
|
|
exp = (s + 1)/b - 1
|
|
fac = fac/(Abs(b) * a1**exp)
|
|
|
|
def tr(l):
|
|
return [a + exp for a in l]
|
|
g1 = meijerg(tr(g1.an), tr(g1.aother), tr(g1.bm), tr(g1.bother), a1*x)
|
|
g2 = meijerg(g2.an, g2.aother, g2.bm, g2.bother, a2*x)
|
|
|
|
from sympy.simplify import powdenest
|
|
return powdenest(fac, polar=True), g1, g2
|
|
|
|
|
|
def _check_antecedents(g1, g2, x):
|
|
""" Return a condition under which the integral theorem applies. """
|
|
# Yes, this is madness.
|
|
# XXX TODO this is a testing *nightmare*
|
|
# NOTE if you update these conditions, please update the documentation as well
|
|
|
|
# The following conditions are found in
|
|
# [P], Section 2.24.1
|
|
#
|
|
# They are also reproduced (verbatim!) at
|
|
# https://functions.wolfram.com/HypergeometricFunctions/MeijerG/21/02/03/
|
|
#
|
|
# Note: k=l=r=alpha=1
|
|
sigma, _ = _get_coeff_exp(g1.argument, x)
|
|
omega, _ = _get_coeff_exp(g2.argument, x)
|
|
s, t, u, v = S([len(g1.bm), len(g1.an), len(g1.ap), len(g1.bq)])
|
|
m, n, p, q = S([len(g2.bm), len(g2.an), len(g2.ap), len(g2.bq)])
|
|
bstar = s + t - (u + v)/2
|
|
cstar = m + n - (p + q)/2
|
|
rho = g1.nu + (u - v)/2 + 1
|
|
mu = g2.nu + (p - q)/2 + 1
|
|
phi = q - p - (v - u)
|
|
eta = 1 - (v - u) - mu - rho
|
|
psi = (pi*(q - m - n) + Abs(unbranched_argument(omega)))/(q - p)
|
|
theta = (pi*(v - s - t) + Abs(unbranched_argument(sigma)))/(v - u)
|
|
|
|
_debug('Checking antecedents:')
|
|
_debugf(' sigma=%s, s=%s, t=%s, u=%s, v=%s, b*=%s, rho=%s',
|
|
(sigma, s, t, u, v, bstar, rho))
|
|
_debugf(' omega=%s, m=%s, n=%s, p=%s, q=%s, c*=%s, mu=%s,',
|
|
(omega, m, n, p, q, cstar, mu))
|
|
_debugf(' phi=%s, eta=%s, psi=%s, theta=%s', (phi, eta, psi, theta))
|
|
|
|
def _c1():
|
|
for g in [g1, g2]:
|
|
for i, j in itertools.product(g.an, g.bm):
|
|
diff = i - j
|
|
if diff.is_integer and diff.is_positive:
|
|
return False
|
|
return True
|
|
c1 = _c1()
|
|
c2 = And(*[re(1 + i + j) > 0 for i in g1.bm for j in g2.bm])
|
|
c3 = And(*[re(1 + i + j) < 1 + 1 for i in g1.an for j in g2.an])
|
|
c4 = And(*[(p - q)*re(1 + i - 1) - re(mu) > Rational(-3, 2) for i in g1.an])
|
|
c5 = And(*[(p - q)*re(1 + i) - re(mu) > Rational(-3, 2) for i in g1.bm])
|
|
c6 = And(*[(u - v)*re(1 + i - 1) - re(rho) > Rational(-3, 2) for i in g2.an])
|
|
c7 = And(*[(u - v)*re(1 + i) - re(rho) > Rational(-3, 2) for i in g2.bm])
|
|
c8 = (Abs(phi) + 2*re((rho - 1)*(q - p) + (v - u)*(q - p) + (mu -
|
|
1)*(v - u)) > 0)
|
|
c9 = (Abs(phi) - 2*re((rho - 1)*(q - p) + (v - u)*(q - p) + (mu -
|
|
1)*(v - u)) > 0)
|
|
c10 = (Abs(unbranched_argument(sigma)) < bstar*pi)
|
|
c11 = Eq(Abs(unbranched_argument(sigma)), bstar*pi)
|
|
c12 = (Abs(unbranched_argument(omega)) < cstar*pi)
|
|
c13 = Eq(Abs(unbranched_argument(omega)), cstar*pi)
|
|
|
|
# The following condition is *not* implemented as stated on the wolfram
|
|
# function site. In the book of Prudnikov there is an additional part
|
|
# (the And involving re()). However, I only have this book in russian, and
|
|
# I don't read any russian. The following condition is what other people
|
|
# have told me it means.
|
|
# Worryingly, it is different from the condition implemented in REDUCE.
|
|
# The REDUCE implementation:
|
|
# https://reduce-algebra.svn.sourceforge.net/svnroot/reduce-algebra/trunk/packages/defint/definta.red
|
|
# (search for tst14)
|
|
# The Wolfram alpha version:
|
|
# https://functions.wolfram.com/HypergeometricFunctions/MeijerG/21/02/03/03/0014/
|
|
z0 = exp(-(bstar + cstar)*pi*S.ImaginaryUnit)
|
|
zos = unpolarify(z0*omega/sigma)
|
|
zso = unpolarify(z0*sigma/omega)
|
|
if zos == 1/zso:
|
|
c14 = And(Eq(phi, 0), bstar + cstar <= 1,
|
|
Or(Ne(zos, 1), re(mu + rho + v - u) < 1,
|
|
re(mu + rho + q - p) < 1))
|
|
else:
|
|
def _cond(z):
|
|
'''Returns True if abs(arg(1-z)) < pi, avoiding arg(0).
|
|
|
|
Explanation
|
|
===========
|
|
|
|
If ``z`` is 1 then arg is NaN. This raises a
|
|
TypeError on `NaN < pi`. Previously this gave `False` so
|
|
this behavior has been hardcoded here but someone should
|
|
check if this NaN is more serious! This NaN is triggered by
|
|
test_meijerint() in test_meijerint.py:
|
|
`meijerint_definite(exp(x), x, 0, I)`
|
|
'''
|
|
return z != 1 and Abs(arg(1 - z)) < pi
|
|
|
|
c14 = And(Eq(phi, 0), bstar - 1 + cstar <= 0,
|
|
Or(And(Ne(zos, 1), _cond(zos)),
|
|
And(re(mu + rho + v - u) < 1, Eq(zos, 1))))
|
|
|
|
c14_alt = And(Eq(phi, 0), cstar - 1 + bstar <= 0,
|
|
Or(And(Ne(zso, 1), _cond(zso)),
|
|
And(re(mu + rho + q - p) < 1, Eq(zso, 1))))
|
|
|
|
# Since r=k=l=1, in our case there is c14_alt which is the same as calling
|
|
# us with (g1, g2) = (g2, g1). The conditions below enumerate all cases
|
|
# (i.e. we don't have to try arguments reversed by hand), and indeed try
|
|
# all symmetric cases. (i.e. whenever there is a condition involving c14,
|
|
# there is also a dual condition which is exactly what we would get when g1,
|
|
# g2 were interchanged, *but c14 was unaltered*).
|
|
# Hence the following seems correct:
|
|
c14 = Or(c14, c14_alt)
|
|
|
|
'''
|
|
When `c15` is NaN (e.g. from `psi` being NaN as happens during
|
|
'test_issue_4992' and/or `theta` is NaN as in 'test_issue_6253',
|
|
both in `test_integrals.py`) the comparison to 0 formerly gave False
|
|
whereas now an error is raised. To keep the old behavior, the value
|
|
of NaN is replaced with False but perhaps a closer look at this condition
|
|
should be made: XXX how should conditions leading to c15=NaN be handled?
|
|
'''
|
|
try:
|
|
lambda_c = (q - p)*Abs(omega)**(1/(q - p))*cos(psi) \
|
|
+ (v - u)*Abs(sigma)**(1/(v - u))*cos(theta)
|
|
# the TypeError might be raised here, e.g. if lambda_c is NaN
|
|
if _eval_cond(lambda_c > 0) != False:
|
|
c15 = (lambda_c > 0)
|
|
else:
|
|
def lambda_s0(c1, c2):
|
|
return c1*(q - p)*Abs(omega)**(1/(q - p))*sin(psi) \
|
|
+ c2*(v - u)*Abs(sigma)**(1/(v - u))*sin(theta)
|
|
lambda_s = Piecewise(
|
|
((lambda_s0(+1, +1)*lambda_s0(-1, -1)),
|
|
And(Eq(unbranched_argument(sigma), 0), Eq(unbranched_argument(omega), 0))),
|
|
(lambda_s0(sign(unbranched_argument(omega)), +1)*lambda_s0(sign(unbranched_argument(omega)), -1),
|
|
And(Eq(unbranched_argument(sigma), 0), Ne(unbranched_argument(omega), 0))),
|
|
(lambda_s0(+1, sign(unbranched_argument(sigma)))*lambda_s0(-1, sign(unbranched_argument(sigma))),
|
|
And(Ne(unbranched_argument(sigma), 0), Eq(unbranched_argument(omega), 0))),
|
|
(lambda_s0(sign(unbranched_argument(omega)), sign(unbranched_argument(sigma))), True))
|
|
tmp = [lambda_c > 0,
|
|
And(Eq(lambda_c, 0), Ne(lambda_s, 0), re(eta) > -1),
|
|
And(Eq(lambda_c, 0), Eq(lambda_s, 0), re(eta) > 0)]
|
|
c15 = Or(*tmp)
|
|
except TypeError:
|
|
c15 = False
|
|
for cond, i in [(c1, 1), (c2, 2), (c3, 3), (c4, 4), (c5, 5), (c6, 6),
|
|
(c7, 7), (c8, 8), (c9, 9), (c10, 10), (c11, 11),
|
|
(c12, 12), (c13, 13), (c14, 14), (c15, 15)]:
|
|
_debugf(' c%s: %s', (i, cond))
|
|
|
|
# We will return Or(*conds)
|
|
conds = []
|
|
|
|
def pr(count):
|
|
_debugf(' case %s: %s', (count, conds[-1]))
|
|
conds += [And(m*n*s*t != 0, bstar.is_positive is True, cstar.is_positive is True, c1, c2, c3, c10,
|
|
c12)] # 1
|
|
pr(1)
|
|
conds += [And(Eq(u, v), Eq(bstar, 0), cstar.is_positive is True, sigma.is_positive is True, re(rho) < 1,
|
|
c1, c2, c3, c12)] # 2
|
|
pr(2)
|
|
conds += [And(Eq(p, q), Eq(cstar, 0), bstar.is_positive is True, omega.is_positive is True, re(mu) < 1,
|
|
c1, c2, c3, c10)] # 3
|
|
pr(3)
|
|
conds += [And(Eq(p, q), Eq(u, v), Eq(bstar, 0), Eq(cstar, 0),
|
|
sigma.is_positive is True, omega.is_positive is True, re(mu) < 1, re(rho) < 1,
|
|
Ne(sigma, omega), c1, c2, c3)] # 4
|
|
pr(4)
|
|
conds += [And(Eq(p, q), Eq(u, v), Eq(bstar, 0), Eq(cstar, 0),
|
|
sigma.is_positive is True, omega.is_positive is True, re(mu + rho) < 1,
|
|
Ne(omega, sigma), c1, c2, c3)] # 5
|
|
pr(5)
|
|
conds += [And(p > q, s.is_positive is True, bstar.is_positive is True, cstar >= 0,
|
|
c1, c2, c3, c5, c10, c13)] # 6
|
|
pr(6)
|
|
conds += [And(p < q, t.is_positive is True, bstar.is_positive is True, cstar >= 0,
|
|
c1, c2, c3, c4, c10, c13)] # 7
|
|
pr(7)
|
|
conds += [And(u > v, m.is_positive is True, cstar.is_positive is True, bstar >= 0,
|
|
c1, c2, c3, c7, c11, c12)] # 8
|
|
pr(8)
|
|
conds += [And(u < v, n.is_positive is True, cstar.is_positive is True, bstar >= 0,
|
|
c1, c2, c3, c6, c11, c12)] # 9
|
|
pr(9)
|
|
conds += [And(p > q, Eq(u, v), Eq(bstar, 0), cstar >= 0, sigma.is_positive is True,
|
|
re(rho) < 1, c1, c2, c3, c5, c13)] # 10
|
|
pr(10)
|
|
conds += [And(p < q, Eq(u, v), Eq(bstar, 0), cstar >= 0, sigma.is_positive is True,
|
|
re(rho) < 1, c1, c2, c3, c4, c13)] # 11
|
|
pr(11)
|
|
conds += [And(Eq(p, q), u > v, bstar >= 0, Eq(cstar, 0), omega.is_positive is True,
|
|
re(mu) < 1, c1, c2, c3, c7, c11)] # 12
|
|
pr(12)
|
|
conds += [And(Eq(p, q), u < v, bstar >= 0, Eq(cstar, 0), omega.is_positive is True,
|
|
re(mu) < 1, c1, c2, c3, c6, c11)] # 13
|
|
pr(13)
|
|
conds += [And(p < q, u > v, bstar >= 0, cstar >= 0,
|
|
c1, c2, c3, c4, c7, c11, c13)] # 14
|
|
pr(14)
|
|
conds += [And(p > q, u < v, bstar >= 0, cstar >= 0,
|
|
c1, c2, c3, c5, c6, c11, c13)] # 15
|
|
pr(15)
|
|
conds += [And(p > q, u > v, bstar >= 0, cstar >= 0,
|
|
c1, c2, c3, c5, c7, c8, c11, c13, c14)] # 16
|
|
pr(16)
|
|
conds += [And(p < q, u < v, bstar >= 0, cstar >= 0,
|
|
c1, c2, c3, c4, c6, c9, c11, c13, c14)] # 17
|
|
pr(17)
|
|
conds += [And(Eq(t, 0), s.is_positive is True, bstar.is_positive is True, phi.is_positive is True, c1, c2, c10)] # 18
|
|
pr(18)
|
|
conds += [And(Eq(s, 0), t.is_positive is True, bstar.is_positive is True, phi.is_negative is True, c1, c3, c10)] # 19
|
|
pr(19)
|
|
conds += [And(Eq(n, 0), m.is_positive is True, cstar.is_positive is True, phi.is_negative is True, c1, c2, c12)] # 20
|
|
pr(20)
|
|
conds += [And(Eq(m, 0), n.is_positive is True, cstar.is_positive is True, phi.is_positive is True, c1, c3, c12)] # 21
|
|
pr(21)
|
|
conds += [And(Eq(s*t, 0), bstar.is_positive is True, cstar.is_positive is True,
|
|
c1, c2, c3, c10, c12)] # 22
|
|
pr(22)
|
|
conds += [And(Eq(m*n, 0), bstar.is_positive is True, cstar.is_positive is True,
|
|
c1, c2, c3, c10, c12)] # 23
|
|
pr(23)
|
|
|
|
# The following case is from [Luke1969]. As far as I can tell, it is *not*
|
|
# covered by Prudnikov's.
|
|
# Let G1 and G2 be the two G-functions. Suppose the integral exists from
|
|
# 0 to a > 0 (this is easy the easy part), that G1 is exponential decay at
|
|
# infinity, and that the mellin transform of G2 exists.
|
|
# Then the integral exists.
|
|
mt1_exists = _check_antecedents_1(g1, x, helper=True)
|
|
mt2_exists = _check_antecedents_1(g2, x, helper=True)
|
|
conds += [And(mt2_exists, Eq(t, 0), u < s, bstar.is_positive is True, c10, c1, c2, c3)]
|
|
pr('E1')
|
|
conds += [And(mt2_exists, Eq(s, 0), v < t, bstar.is_positive is True, c10, c1, c2, c3)]
|
|
pr('E2')
|
|
conds += [And(mt1_exists, Eq(n, 0), p < m, cstar.is_positive is True, c12, c1, c2, c3)]
|
|
pr('E3')
|
|
conds += [And(mt1_exists, Eq(m, 0), q < n, cstar.is_positive is True, c12, c1, c2, c3)]
|
|
pr('E4')
|
|
|
|
# Let's short-circuit if this worked ...
|
|
# the rest is corner-cases and terrible to read.
|
|
r = Or(*conds)
|
|
if _eval_cond(r) != False:
|
|
return r
|
|
|
|
conds += [And(m + n > p, Eq(t, 0), Eq(phi, 0), s.is_positive is True, bstar.is_positive is True, cstar.is_negative is True,
|
|
Abs(unbranched_argument(omega)) < (m + n - p + 1)*pi,
|
|
c1, c2, c10, c14, c15)] # 24
|
|
pr(24)
|
|
conds += [And(m + n > q, Eq(s, 0), Eq(phi, 0), t.is_positive is True, bstar.is_positive is True, cstar.is_negative is True,
|
|
Abs(unbranched_argument(omega)) < (m + n - q + 1)*pi,
|
|
c1, c3, c10, c14, c15)] # 25
|
|
pr(25)
|
|
conds += [And(Eq(p, q - 1), Eq(t, 0), Eq(phi, 0), s.is_positive is True, bstar.is_positive is True,
|
|
cstar >= 0, cstar*pi < Abs(unbranched_argument(omega)),
|
|
c1, c2, c10, c14, c15)] # 26
|
|
pr(26)
|
|
conds += [And(Eq(p, q + 1), Eq(s, 0), Eq(phi, 0), t.is_positive is True, bstar.is_positive is True,
|
|
cstar >= 0, cstar*pi < Abs(unbranched_argument(omega)),
|
|
c1, c3, c10, c14, c15)] # 27
|
|
pr(27)
|
|
conds += [And(p < q - 1, Eq(t, 0), Eq(phi, 0), s.is_positive is True, bstar.is_positive is True,
|
|
cstar >= 0, cstar*pi < Abs(unbranched_argument(omega)),
|
|
Abs(unbranched_argument(omega)) < (m + n - p + 1)*pi,
|
|
c1, c2, c10, c14, c15)] # 28
|
|
pr(28)
|
|
conds += [And(
|
|
p > q + 1, Eq(s, 0), Eq(phi, 0), t.is_positive is True, bstar.is_positive is True, cstar >= 0,
|
|
cstar*pi < Abs(unbranched_argument(omega)),
|
|
Abs(unbranched_argument(omega)) < (m + n - q + 1)*pi,
|
|
c1, c3, c10, c14, c15)] # 29
|
|
pr(29)
|
|
conds += [And(Eq(n, 0), Eq(phi, 0), s + t > 0, m.is_positive is True, cstar.is_positive is True, bstar.is_negative is True,
|
|
Abs(unbranched_argument(sigma)) < (s + t - u + 1)*pi,
|
|
c1, c2, c12, c14, c15)] # 30
|
|
pr(30)
|
|
conds += [And(Eq(m, 0), Eq(phi, 0), s + t > v, n.is_positive is True, cstar.is_positive is True, bstar.is_negative is True,
|
|
Abs(unbranched_argument(sigma)) < (s + t - v + 1)*pi,
|
|
c1, c3, c12, c14, c15)] # 31
|
|
pr(31)
|
|
conds += [And(Eq(n, 0), Eq(phi, 0), Eq(u, v - 1), m.is_positive is True, cstar.is_positive is True,
|
|
bstar >= 0, bstar*pi < Abs(unbranched_argument(sigma)),
|
|
Abs(unbranched_argument(sigma)) < (bstar + 1)*pi,
|
|
c1, c2, c12, c14, c15)] # 32
|
|
pr(32)
|
|
conds += [And(Eq(m, 0), Eq(phi, 0), Eq(u, v + 1), n.is_positive is True, cstar.is_positive is True,
|
|
bstar >= 0, bstar*pi < Abs(unbranched_argument(sigma)),
|
|
Abs(unbranched_argument(sigma)) < (bstar + 1)*pi,
|
|
c1, c3, c12, c14, c15)] # 33
|
|
pr(33)
|
|
conds += [And(
|
|
Eq(n, 0), Eq(phi, 0), u < v - 1, m.is_positive is True, cstar.is_positive is True, bstar >= 0,
|
|
bstar*pi < Abs(unbranched_argument(sigma)),
|
|
Abs(unbranched_argument(sigma)) < (s + t - u + 1)*pi,
|
|
c1, c2, c12, c14, c15)] # 34
|
|
pr(34)
|
|
conds += [And(
|
|
Eq(m, 0), Eq(phi, 0), u > v + 1, n.is_positive is True, cstar.is_positive is True, bstar >= 0,
|
|
bstar*pi < Abs(unbranched_argument(sigma)),
|
|
Abs(unbranched_argument(sigma)) < (s + t - v + 1)*pi,
|
|
c1, c3, c12, c14, c15)] # 35
|
|
pr(35)
|
|
|
|
return Or(*conds)
|
|
|
|
# NOTE An alternative, but as far as I can tell weaker, set of conditions
|
|
# can be found in [L, section 5.6.2].
|
|
|
|
|
|
def _int0oo(g1, g2, x):
|
|
"""
|
|
Express integral from zero to infinity g1*g2 using a G function,
|
|
assuming the necessary conditions are fulfilled.
|
|
|
|
Examples
|
|
========
|
|
|
|
>>> from sympy.integrals.meijerint import _int0oo
|
|
>>> from sympy.abc import s, t, m
|
|
>>> from sympy import meijerg, S
|
|
>>> g1 = meijerg([], [], [-S(1)/2, 0], [], s**2*t/4)
|
|
>>> g2 = meijerg([], [], [m/2], [-m/2], t/4)
|
|
>>> _int0oo(g1, g2, t)
|
|
4*meijerg(((1/2, 0), ()), ((m/2,), (-m/2,)), s**(-2))/s**2
|
|
"""
|
|
# See: [L, section 5.6.2, equation (1)]
|
|
eta, _ = _get_coeff_exp(g1.argument, x)
|
|
omega, _ = _get_coeff_exp(g2.argument, x)
|
|
|
|
def neg(l):
|
|
return [-x for x in l]
|
|
a1 = neg(g1.bm) + list(g2.an)
|
|
a2 = list(g2.aother) + neg(g1.bother)
|
|
b1 = neg(g1.an) + list(g2.bm)
|
|
b2 = list(g2.bother) + neg(g1.aother)
|
|
return meijerg(a1, a2, b1, b2, omega/eta)/eta
|
|
|
|
|
|
def _rewrite_inversion(fac, po, g, x):
|
|
""" Absorb ``po`` == x**s into g. """
|
|
_, s = _get_coeff_exp(po, x)
|
|
a, b = _get_coeff_exp(g.argument, x)
|
|
|
|
def tr(l):
|
|
return [t + s/b for t in l]
|
|
from sympy.simplify import powdenest
|
|
return (powdenest(fac/a**(s/b), polar=True),
|
|
meijerg(tr(g.an), tr(g.aother), tr(g.bm), tr(g.bother), g.argument))
|
|
|
|
|
|
def _check_antecedents_inversion(g, x):
|
|
""" Check antecedents for the laplace inversion integral. """
|
|
_debug('Checking antecedents for inversion:')
|
|
z = g.argument
|
|
_, e = _get_coeff_exp(z, x)
|
|
if e < 0:
|
|
_debug(' Flipping G.')
|
|
# We want to assume that argument gets large as |x| -> oo
|
|
return _check_antecedents_inversion(_flip_g(g), x)
|
|
|
|
def statement_half(a, b, c, z, plus):
|
|
coeff, exponent = _get_coeff_exp(z, x)
|
|
a *= exponent
|
|
b *= coeff**c
|
|
c *= exponent
|
|
conds = []
|
|
wp = b*exp(S.ImaginaryUnit*re(c)*pi/2)
|
|
wm = b*exp(-S.ImaginaryUnit*re(c)*pi/2)
|
|
if plus:
|
|
w = wp
|
|
else:
|
|
w = wm
|
|
conds += [And(Or(Eq(b, 0), re(c) <= 0), re(a) <= -1)]
|
|
conds += [And(Ne(b, 0), Eq(im(c), 0), re(c) > 0, re(w) < 0)]
|
|
conds += [And(Ne(b, 0), Eq(im(c), 0), re(c) > 0, re(w) <= 0,
|
|
re(a) <= -1)]
|
|
return Or(*conds)
|
|
|
|
def statement(a, b, c, z):
|
|
""" Provide a convergence statement for z**a * exp(b*z**c),
|
|
c/f sphinx docs. """
|
|
return And(statement_half(a, b, c, z, True),
|
|
statement_half(a, b, c, z, False))
|
|
|
|
# Notations from [L], section 5.7-10
|
|
m, n, p, q = S([len(g.bm), len(g.an), len(g.ap), len(g.bq)])
|
|
tau = m + n - p
|
|
nu = q - m - n
|
|
rho = (tau - nu)/2
|
|
sigma = q - p
|
|
if sigma == 1:
|
|
epsilon = S.Half
|
|
elif sigma > 1:
|
|
epsilon = 1
|
|
else:
|
|
epsilon = S.NaN
|
|
theta = ((1 - sigma)/2 + Add(*g.bq) - Add(*g.ap))/sigma
|
|
delta = g.delta
|
|
_debugf(' m=%s, n=%s, p=%s, q=%s, tau=%s, nu=%s, rho=%s, sigma=%s',
|
|
(m, n, p, q, tau, nu, rho, sigma))
|
|
_debugf(' epsilon=%s, theta=%s, delta=%s', (epsilon, theta, delta))
|
|
|
|
# First check if the computation is valid.
|
|
if not (g.delta >= e/2 or (p >= 1 and p >= q)):
|
|
_debug(' Computation not valid for these parameters.')
|
|
return False
|
|
|
|
# Now check if the inversion integral exists.
|
|
|
|
# Test "condition A"
|
|
for a, b in itertools.product(g.an, g.bm):
|
|
if (a - b).is_integer and a > b:
|
|
_debug(' Not a valid G function.')
|
|
return False
|
|
|
|
# There are two cases. If p >= q, we can directly use a slater expansion
|
|
# like [L], 5.2 (11). Note in particular that the asymptotics of such an
|
|
# expansion even hold when some of the parameters differ by integers, i.e.
|
|
# the formula itself would not be valid! (b/c G functions are cts. in their
|
|
# parameters)
|
|
# When p < q, we need to use the theorems of [L], 5.10.
|
|
|
|
if p >= q:
|
|
_debug(' Using asymptotic Slater expansion.')
|
|
return And(*[statement(a - 1, 0, 0, z) for a in g.an])
|
|
|
|
def E(z):
|
|
return And(*[statement(a - 1, 0, 0, z) for a in g.an])
|
|
|
|
def H(z):
|
|
return statement(theta, -sigma, 1/sigma, z)
|
|
|
|
def Hp(z):
|
|
return statement_half(theta, -sigma, 1/sigma, z, True)
|
|
|
|
def Hm(z):
|
|
return statement_half(theta, -sigma, 1/sigma, z, False)
|
|
|
|
# [L], section 5.10
|
|
conds = []
|
|
# Theorem 1 -- p < q from test above
|
|
conds += [And(1 <= n, 1 <= m, rho*pi - delta >= pi/2, delta > 0,
|
|
E(z*exp(S.ImaginaryUnit*pi*(nu + 1))))]
|
|
# Theorem 2, statements (2) and (3)
|
|
conds += [And(p + 1 <= m, m + 1 <= q, delta > 0, delta < pi/2, n == 0,
|
|
(m - p + 1)*pi - delta >= pi/2,
|
|
Hp(z*exp(S.ImaginaryUnit*pi*(q - m))),
|
|
Hm(z*exp(-S.ImaginaryUnit*pi*(q - m))))]
|
|
# Theorem 2, statement (5) -- p < q from test above
|
|
conds += [And(m == q, n == 0, delta > 0,
|
|
(sigma + epsilon)*pi - delta >= pi/2, H(z))]
|
|
# Theorem 3, statements (6) and (7)
|
|
conds += [And(Or(And(p <= q - 2, 1 <= tau, tau <= sigma/2),
|
|
And(p + 1 <= m + n, m + n <= (p + q)/2)),
|
|
delta > 0, delta < pi/2, (tau + 1)*pi - delta >= pi/2,
|
|
Hp(z*exp(S.ImaginaryUnit*pi*nu)),
|
|
Hm(z*exp(-S.ImaginaryUnit*pi*nu)))]
|
|
# Theorem 4, statements (10) and (11) -- p < q from test above
|
|
conds += [And(1 <= m, rho > 0, delta > 0, delta + rho*pi < pi/2,
|
|
(tau + epsilon)*pi - delta >= pi/2,
|
|
Hp(z*exp(S.ImaginaryUnit*pi*nu)),
|
|
Hm(z*exp(-S.ImaginaryUnit*pi*nu)))]
|
|
# Trivial case
|
|
conds += [m == 0]
|
|
|
|
# TODO
|
|
# Theorem 5 is quite general
|
|
# Theorem 6 contains special cases for q=p+1
|
|
|
|
return Or(*conds)
|
|
|
|
|
|
def _int_inversion(g, x, t):
|
|
"""
|
|
Compute the laplace inversion integral, assuming the formula applies.
|
|
"""
|
|
b, a = _get_coeff_exp(g.argument, x)
|
|
C, g = _inflate_fox_h(meijerg(g.an, g.aother, g.bm, g.bother, b/t**a), -a)
|
|
return C/t*g
|
|
|
|
|
|
####################################################################
|
|
# Finally, the real meat.
|
|
####################################################################
|
|
|
|
_lookup_table = None
|
|
|
|
|
|
@cacheit
|
|
@timeit
|
|
def _rewrite_single(f, x, recursive=True):
|
|
"""
|
|
Try to rewrite f as a sum of single G functions of the form
|
|
C*x**s*G(a*x**b), where b is a rational number and C is independent of x.
|
|
We guarantee that result.argument.as_coeff_mul(x) returns (a, (x**b,))
|
|
or (a, ()).
|
|
Returns a list of tuples (C, s, G) and a condition cond.
|
|
Returns None on failure.
|
|
"""
|
|
from .transforms import (mellin_transform, inverse_mellin_transform,
|
|
IntegralTransformError, MellinTransformStripError)
|
|
|
|
global _lookup_table
|
|
if not _lookup_table:
|
|
_lookup_table = {}
|
|
_create_lookup_table(_lookup_table)
|
|
|
|
if isinstance(f, meijerg):
|
|
coeff, m = factor(f.argument, x).as_coeff_mul(x)
|
|
if len(m) > 1:
|
|
return None
|
|
m = m[0]
|
|
if m.is_Pow:
|
|
if m.base != x or not m.exp.is_Rational:
|
|
return None
|
|
elif m != x:
|
|
return None
|
|
return [(1, 0, meijerg(f.an, f.aother, f.bm, f.bother, coeff*m))], True
|
|
|
|
f_ = f
|
|
f = f.subs(x, z)
|
|
t = _mytype(f, z)
|
|
if t in _lookup_table:
|
|
l = _lookup_table[t]
|
|
for formula, terms, cond, hint in l:
|
|
subs = f.match(formula, old=True)
|
|
if subs:
|
|
subs_ = {}
|
|
for fro, to in subs.items():
|
|
subs_[fro] = unpolarify(polarify(to, lift=True),
|
|
exponents_only=True)
|
|
subs = subs_
|
|
if not isinstance(hint, bool):
|
|
hint = hint.subs(subs)
|
|
if hint == False:
|
|
continue
|
|
if not isinstance(cond, (bool, BooleanAtom)):
|
|
cond = unpolarify(cond.subs(subs))
|
|
if _eval_cond(cond) == False:
|
|
continue
|
|
if not isinstance(terms, list):
|
|
terms = terms(subs)
|
|
res = []
|
|
for fac, g in terms:
|
|
r1 = _get_coeff_exp(unpolarify(fac.subs(subs).subs(z, x),
|
|
exponents_only=True), x)
|
|
try:
|
|
g = g.subs(subs).subs(z, x)
|
|
except ValueError:
|
|
continue
|
|
# NOTE these substitutions can in principle introduce oo,
|
|
# zoo and other absurdities. It shouldn't matter,
|
|
# but better be safe.
|
|
if Tuple(*(r1 + (g,))).has(S.Infinity, S.ComplexInfinity, S.NegativeInfinity):
|
|
continue
|
|
g = meijerg(g.an, g.aother, g.bm, g.bother,
|
|
unpolarify(g.argument, exponents_only=True))
|
|
res.append(r1 + (g,))
|
|
if res:
|
|
return res, cond
|
|
|
|
# try recursive mellin transform
|
|
if not recursive:
|
|
return None
|
|
_debug('Trying recursive Mellin transform method.')
|
|
|
|
def my_imt(F, s, x, strip):
|
|
""" Calling simplify() all the time is slow and not helpful, since
|
|
most of the time it only factors things in a way that has to be
|
|
un-done anyway. But sometimes it can remove apparent poles. """
|
|
# XXX should this be in inverse_mellin_transform?
|
|
try:
|
|
return inverse_mellin_transform(F, s, x, strip,
|
|
as_meijerg=True, needeval=True)
|
|
except MellinTransformStripError:
|
|
from sympy.simplify import simplify
|
|
return inverse_mellin_transform(
|
|
simplify(cancel(expand(F))), s, x, strip,
|
|
as_meijerg=True, needeval=True)
|
|
f = f_
|
|
s = _dummy('s', 'rewrite-single', f)
|
|
# to avoid infinite recursion, we have to force the two g functions case
|
|
|
|
def my_integrator(f, x):
|
|
r = _meijerint_definite_4(f, x, only_double=True)
|
|
if r is not None:
|
|
from sympy.simplify import hyperexpand
|
|
res, cond = r
|
|
res = _my_unpolarify(hyperexpand(res, rewrite='nonrepsmall'))
|
|
return Piecewise((res, cond),
|
|
(Integral(f, (x, S.Zero, S.Infinity)), True))
|
|
return Integral(f, (x, S.Zero, S.Infinity))
|
|
try:
|
|
F, strip, _ = mellin_transform(f, x, s, integrator=my_integrator,
|
|
simplify=False, needeval=True)
|
|
g = my_imt(F, s, x, strip)
|
|
except IntegralTransformError:
|
|
g = None
|
|
if g is None:
|
|
# We try to find an expression by analytic continuation.
|
|
# (also if the dummy is already in the expression, there is no point in
|
|
# putting in another one)
|
|
a = _dummy_('a', 'rewrite-single')
|
|
if a not in f.free_symbols and _is_analytic(f, x):
|
|
try:
|
|
F, strip, _ = mellin_transform(f.subs(x, a*x), x, s,
|
|
integrator=my_integrator,
|
|
needeval=True, simplify=False)
|
|
g = my_imt(F, s, x, strip).subs(a, 1)
|
|
except IntegralTransformError:
|
|
g = None
|
|
if g is None or g.has(S.Infinity, S.NaN, S.ComplexInfinity):
|
|
_debug('Recursive Mellin transform failed.')
|
|
return None
|
|
args = Add.make_args(g)
|
|
res = []
|
|
for f in args:
|
|
c, m = f.as_coeff_mul(x)
|
|
if len(m) > 1:
|
|
raise NotImplementedError('Unexpected form...')
|
|
g = m[0]
|
|
a, b = _get_coeff_exp(g.argument, x)
|
|
res += [(c, 0, meijerg(g.an, g.aother, g.bm, g.bother,
|
|
unpolarify(polarify(
|
|
a, lift=True), exponents_only=True)
|
|
*x**b))]
|
|
_debug('Recursive Mellin transform worked:', g)
|
|
return res, True
|
|
|
|
|
|
def _rewrite1(f, x, recursive=True):
|
|
"""
|
|
Try to rewrite ``f`` using a (sum of) single G functions with argument a*x**b.
|
|
Return fac, po, g such that f = fac*po*g, fac is independent of ``x``.
|
|
and po = x**s.
|
|
Here g is a result from _rewrite_single.
|
|
Return None on failure.
|
|
"""
|
|
fac, po, g = _split_mul(f, x)
|
|
g = _rewrite_single(g, x, recursive)
|
|
if g:
|
|
return fac, po, g[0], g[1]
|
|
|
|
|
|
def _rewrite2(f, x):
|
|
"""
|
|
Try to rewrite ``f`` as a product of two G functions of arguments a*x**b.
|
|
Return fac, po, g1, g2 such that f = fac*po*g1*g2, where fac is
|
|
independent of x and po is x**s.
|
|
Here g1 and g2 are results of _rewrite_single.
|
|
Returns None on failure.
|
|
"""
|
|
fac, po, g = _split_mul(f, x)
|
|
if any(_rewrite_single(expr, x, False) is None for expr in _mul_args(g)):
|
|
return None
|
|
l = _mul_as_two_parts(g)
|
|
if not l:
|
|
return None
|
|
l = list(ordered(l, [
|
|
lambda p: max(len(_exponents(p[0], x)), len(_exponents(p[1], x))),
|
|
lambda p: max(len(_functions(p[0], x)), len(_functions(p[1], x))),
|
|
lambda p: max(len(_find_splitting_points(p[0], x)),
|
|
len(_find_splitting_points(p[1], x)))]))
|
|
|
|
for recursive, (fac1, fac2) in itertools.product((False, True), l):
|
|
g1 = _rewrite_single(fac1, x, recursive)
|
|
g2 = _rewrite_single(fac2, x, recursive)
|
|
if g1 and g2:
|
|
cond = And(g1[1], g2[1])
|
|
if cond != False:
|
|
return fac, po, g1[0], g2[0], cond
|
|
|
|
|
|
def meijerint_indefinite(f, x):
|
|
"""
|
|
Compute an indefinite integral of ``f`` by rewriting it as a G function.
|
|
|
|
Examples
|
|
========
|
|
|
|
>>> from sympy.integrals.meijerint import meijerint_indefinite
|
|
>>> from sympy import sin
|
|
>>> from sympy.abc import x
|
|
>>> meijerint_indefinite(sin(x), x)
|
|
-cos(x)
|
|
"""
|
|
f = sympify(f)
|
|
results = []
|
|
for a in sorted(_find_splitting_points(f, x) | {S.Zero}, key=default_sort_key):
|
|
res = _meijerint_indefinite_1(f.subs(x, x + a), x)
|
|
if not res:
|
|
continue
|
|
res = res.subs(x, x - a)
|
|
if _has(res, hyper, meijerg):
|
|
results.append(res)
|
|
else:
|
|
return res
|
|
if f.has(HyperbolicFunction):
|
|
_debug('Try rewriting hyperbolics in terms of exp.')
|
|
rv = meijerint_indefinite(
|
|
_rewrite_hyperbolics_as_exp(f), x)
|
|
if rv:
|
|
if not isinstance(rv, list):
|
|
from sympy.simplify.radsimp import collect
|
|
return collect(factor_terms(rv), rv.atoms(exp))
|
|
results.extend(rv)
|
|
if results:
|
|
return next(ordered(results))
|
|
|
|
|
|
def _meijerint_indefinite_1(f, x):
|
|
""" Helper that does not attempt any substitution. """
|
|
_debug('Trying to compute the indefinite integral of', f, 'wrt', x)
|
|
from sympy.simplify import hyperexpand, powdenest
|
|
|
|
gs = _rewrite1(f, x)
|
|
if gs is None:
|
|
# Note: the code that calls us will do expand() and try again
|
|
return None
|
|
|
|
fac, po, gl, cond = gs
|
|
_debug(' could rewrite:', gs)
|
|
res = S.Zero
|
|
for C, s, g in gl:
|
|
a, b = _get_coeff_exp(g.argument, x)
|
|
_, c = _get_coeff_exp(po, x)
|
|
c += s
|
|
|
|
# we do a substitution t=a*x**b, get integrand fac*t**rho*g
|
|
fac_ = fac * C / (b*a**((1 + c)/b))
|
|
rho = (c + 1)/b - 1
|
|
|
|
# we now use t**rho*G(params, t) = G(params + rho, t)
|
|
# [L, page 150, equation (4)]
|
|
# and integral G(params, t) dt = G(1, params+1, 0, t)
|
|
# (or a similar expression with 1 and 0 exchanged ... pick the one
|
|
# which yields a well-defined function)
|
|
# [R, section 5]
|
|
# (Note that this dummy will immediately go away again, so we
|
|
# can safely pass S.One for ``expr``.)
|
|
t = _dummy('t', 'meijerint-indefinite', S.One)
|
|
|
|
def tr(p):
|
|
return [a + rho + 1 for a in p]
|
|
if any(b.is_integer and (b <= 0) == True for b in tr(g.bm)):
|
|
r = -meijerg(
|
|
tr(g.an), tr(g.aother) + [1], tr(g.bm) + [0], tr(g.bother), t)
|
|
else:
|
|
r = meijerg(
|
|
tr(g.an) + [1], tr(g.aother), tr(g.bm), tr(g.bother) + [0], t)
|
|
# The antiderivative is most often expected to be defined
|
|
# in the neighborhood of x = 0.
|
|
if b.is_extended_nonnegative and not f.subs(x, 0).has(S.NaN, S.ComplexInfinity):
|
|
place = 0 # Assume we can expand at zero
|
|
else:
|
|
place = None
|
|
r = hyperexpand(r.subs(t, a*x**b), place=place)
|
|
|
|
# now substitute back
|
|
# Note: we really do want the powers of x to combine.
|
|
res += powdenest(fac_*r, polar=True)
|
|
|
|
def _clean(res):
|
|
"""This multiplies out superfluous powers of x we created, and chops off
|
|
constants:
|
|
|
|
>> _clean(x*(exp(x)/x - 1/x) + 3)
|
|
exp(x)
|
|
|
|
cancel is used before mul_expand since it is possible for an
|
|
expression to have an additive constant that does not become isolated
|
|
with simple expansion. Such a situation was identified in issue 6369:
|
|
|
|
Examples
|
|
========
|
|
|
|
>>> from sympy import sqrt, cancel
|
|
>>> from sympy.abc import x
|
|
>>> a = sqrt(2*x + 1)
|
|
>>> bad = (3*x*a**5 + 2*x - a**5 + 1)/a**2
|
|
>>> bad.expand().as_independent(x)[0]
|
|
0
|
|
>>> cancel(bad).expand().as_independent(x)[0]
|
|
1
|
|
"""
|
|
res = expand_mul(cancel(res), deep=False)
|
|
return Add._from_args(res.as_coeff_add(x)[1])
|
|
|
|
res = piecewise_fold(res, evaluate=None)
|
|
if res.is_Piecewise:
|
|
newargs = []
|
|
for e, c in res.args:
|
|
e = _my_unpolarify(_clean(e))
|
|
newargs += [(e, c)]
|
|
res = Piecewise(*newargs, evaluate=False)
|
|
else:
|
|
res = _my_unpolarify(_clean(res))
|
|
return Piecewise((res, _my_unpolarify(cond)), (Integral(f, x), True))
|
|
|
|
|
|
@timeit
|
|
def meijerint_definite(f, x, a, b):
|
|
"""
|
|
Integrate ``f`` over the interval [``a``, ``b``], by rewriting it as a product
|
|
of two G functions, or as a single G function.
|
|
|
|
Return res, cond, where cond are convergence conditions.
|
|
|
|
Examples
|
|
========
|
|
|
|
>>> from sympy.integrals.meijerint import meijerint_definite
|
|
>>> from sympy import exp, oo
|
|
>>> from sympy.abc import x
|
|
>>> meijerint_definite(exp(-x**2), x, -oo, oo)
|
|
(sqrt(pi), True)
|
|
|
|
This function is implemented as a succession of functions
|
|
meijerint_definite, _meijerint_definite_2, _meijerint_definite_3,
|
|
_meijerint_definite_4. Each function in the list calls the next one
|
|
(presumably) several times. This means that calling meijerint_definite
|
|
can be very costly.
|
|
"""
|
|
# This consists of three steps:
|
|
# 1) Change the integration limits to 0, oo
|
|
# 2) Rewrite in terms of G functions
|
|
# 3) Evaluate the integral
|
|
#
|
|
# There are usually several ways of doing this, and we want to try all.
|
|
# This function does (1), calls _meijerint_definite_2 for step (2).
|
|
_debugf('Integrating %s wrt %s from %s to %s.', (f, x, a, b))
|
|
f = sympify(f)
|
|
if f.has(DiracDelta):
|
|
_debug('Integrand has DiracDelta terms - giving up.')
|
|
return None
|
|
|
|
if f.has(SingularityFunction):
|
|
_debug('Integrand has Singularity Function terms - giving up.')
|
|
return None
|
|
|
|
f_, x_, a_, b_ = f, x, a, b
|
|
|
|
# Let's use a dummy in case any of the boundaries has x.
|
|
d = Dummy('x')
|
|
f = f.subs(x, d)
|
|
x = d
|
|
|
|
if a == b:
|
|
return (S.Zero, True)
|
|
|
|
results = []
|
|
if a is S.NegativeInfinity and b is not S.Infinity:
|
|
return meijerint_definite(f.subs(x, -x), x, -b, -a)
|
|
|
|
elif a is S.NegativeInfinity:
|
|
# Integrating -oo to oo. We need to find a place to split the integral.
|
|
_debug(' Integrating -oo to +oo.')
|
|
innermost = _find_splitting_points(f, x)
|
|
_debug(' Sensible splitting points:', innermost)
|
|
for c in sorted(innermost, key=default_sort_key, reverse=True) + [S.Zero]:
|
|
_debug(' Trying to split at', c)
|
|
if not c.is_extended_real:
|
|
_debug(' Non-real splitting point.')
|
|
continue
|
|
res1 = _meijerint_definite_2(f.subs(x, x + c), x)
|
|
if res1 is None:
|
|
_debug(' But could not compute first integral.')
|
|
continue
|
|
res2 = _meijerint_definite_2(f.subs(x, c - x), x)
|
|
if res2 is None:
|
|
_debug(' But could not compute second integral.')
|
|
continue
|
|
res1, cond1 = res1
|
|
res2, cond2 = res2
|
|
cond = _condsimp(And(cond1, cond2))
|
|
if cond == False:
|
|
_debug(' But combined condition is always false.')
|
|
continue
|
|
res = res1 + res2
|
|
return res, cond
|
|
|
|
elif a is S.Infinity:
|
|
res = meijerint_definite(f, x, b, S.Infinity)
|
|
return -res[0], res[1]
|
|
|
|
elif (a, b) == (S.Zero, S.Infinity):
|
|
# This is a common case - try it directly first.
|
|
res = _meijerint_definite_2(f, x)
|
|
if res:
|
|
if _has(res[0], meijerg):
|
|
results.append(res)
|
|
else:
|
|
return res
|
|
|
|
else:
|
|
if b is S.Infinity:
|
|
for split in _find_splitting_points(f, x):
|
|
if (a - split >= 0) == True:
|
|
_debugf('Trying x -> x + %s', split)
|
|
res = _meijerint_definite_2(f.subs(x, x + split)
|
|
*Heaviside(x + split - a), x)
|
|
if res:
|
|
if _has(res[0], meijerg):
|
|
results.append(res)
|
|
else:
|
|
return res
|
|
|
|
f = f.subs(x, x + a)
|
|
b = b - a
|
|
a = 0
|
|
if b is not S.Infinity:
|
|
phi = exp(S.ImaginaryUnit*arg(b))
|
|
b = Abs(b)
|
|
f = f.subs(x, phi*x)
|
|
f *= Heaviside(b - x)*phi
|
|
b = S.Infinity
|
|
|
|
_debug('Changed limits to', a, b)
|
|
_debug('Changed function to', f)
|
|
res = _meijerint_definite_2(f, x)
|
|
if res:
|
|
if _has(res[0], meijerg):
|
|
results.append(res)
|
|
else:
|
|
return res
|
|
if f_.has(HyperbolicFunction):
|
|
_debug('Try rewriting hyperbolics in terms of exp.')
|
|
rv = meijerint_definite(
|
|
_rewrite_hyperbolics_as_exp(f_), x_, a_, b_)
|
|
if rv:
|
|
if not isinstance(rv, list):
|
|
from sympy.simplify.radsimp import collect
|
|
rv = (collect(factor_terms(rv[0]), rv[0].atoms(exp)),) + rv[1:]
|
|
return rv
|
|
results.extend(rv)
|
|
if results:
|
|
return next(ordered(results))
|
|
|
|
|
|
def _guess_expansion(f, x):
|
|
""" Try to guess sensible rewritings for integrand f(x). """
|
|
res = [(f, 'original integrand')]
|
|
|
|
orig = res[-1][0]
|
|
saw = {orig}
|
|
expanded = expand_mul(orig)
|
|
if expanded not in saw:
|
|
res += [(expanded, 'expand_mul')]
|
|
saw.add(expanded)
|
|
|
|
expanded = expand(orig)
|
|
if expanded not in saw:
|
|
res += [(expanded, 'expand')]
|
|
saw.add(expanded)
|
|
|
|
if orig.has(TrigonometricFunction, HyperbolicFunction):
|
|
expanded = expand_mul(expand_trig(orig))
|
|
if expanded not in saw:
|
|
res += [(expanded, 'expand_trig, expand_mul')]
|
|
saw.add(expanded)
|
|
|
|
if orig.has(cos, sin):
|
|
from sympy.simplify.fu import sincos_to_sum
|
|
reduced = sincos_to_sum(orig)
|
|
if reduced not in saw:
|
|
res += [(reduced, 'trig power reduction')]
|
|
saw.add(reduced)
|
|
|
|
return res
|
|
|
|
|
|
def _meijerint_definite_2(f, x):
|
|
"""
|
|
Try to integrate f dx from zero to infinity.
|
|
|
|
The body of this function computes various 'simplifications'
|
|
f1, f2, ... of f (e.g. by calling expand_mul(), trigexpand()
|
|
- see _guess_expansion) and calls _meijerint_definite_3 with each of
|
|
these in succession.
|
|
If _meijerint_definite_3 succeeds with any of the simplified functions,
|
|
returns this result.
|
|
"""
|
|
# This function does preparation for (2), calls
|
|
# _meijerint_definite_3 for (2) and (3) combined.
|
|
|
|
# use a positive dummy - we integrate from 0 to oo
|
|
# XXX if a nonnegative symbol is used there will be test failures
|
|
dummy = _dummy('x', 'meijerint-definite2', f, positive=True)
|
|
f = f.subs(x, dummy)
|
|
x = dummy
|
|
|
|
if f == 0:
|
|
return S.Zero, True
|
|
|
|
for g, explanation in _guess_expansion(f, x):
|
|
_debug('Trying', explanation)
|
|
res = _meijerint_definite_3(g, x)
|
|
if res:
|
|
return res
|
|
|
|
|
|
def _meijerint_definite_3(f, x):
|
|
"""
|
|
Try to integrate f dx from zero to infinity.
|
|
|
|
This function calls _meijerint_definite_4 to try to compute the
|
|
integral. If this fails, it tries using linearity.
|
|
"""
|
|
res = _meijerint_definite_4(f, x)
|
|
if res and res[1] != False:
|
|
return res
|
|
if f.is_Add:
|
|
_debug('Expanding and evaluating all terms.')
|
|
ress = [_meijerint_definite_4(g, x) for g in f.args]
|
|
if all(r is not None for r in ress):
|
|
conds = []
|
|
res = S.Zero
|
|
for r, c in ress:
|
|
res += r
|
|
conds += [c]
|
|
c = And(*conds)
|
|
if c != False:
|
|
return res, c
|
|
|
|
|
|
def _my_unpolarify(f):
|
|
return _eval_cond(unpolarify(f))
|
|
|
|
|
|
@timeit
|
|
def _meijerint_definite_4(f, x, only_double=False):
|
|
"""
|
|
Try to integrate f dx from zero to infinity.
|
|
|
|
Explanation
|
|
===========
|
|
|
|
This function tries to apply the integration theorems found in literature,
|
|
i.e. it tries to rewrite f as either one or a product of two G-functions.
|
|
|
|
The parameter ``only_double`` is used internally in the recursive algorithm
|
|
to disable trying to rewrite f as a single G-function.
|
|
"""
|
|
from sympy.simplify import hyperexpand
|
|
# This function does (2) and (3)
|
|
_debug('Integrating', f)
|
|
# Try single G function.
|
|
if not only_double:
|
|
gs = _rewrite1(f, x, recursive=False)
|
|
if gs is not None:
|
|
fac, po, g, cond = gs
|
|
_debug('Could rewrite as single G function:', fac, po, g)
|
|
res = S.Zero
|
|
for C, s, f in g:
|
|
if C == 0:
|
|
continue
|
|
C, f = _rewrite_saxena_1(fac*C, po*x**s, f, x)
|
|
res += C*_int0oo_1(f, x)
|
|
cond = And(cond, _check_antecedents_1(f, x))
|
|
if cond == False:
|
|
break
|
|
cond = _my_unpolarify(cond)
|
|
if cond == False:
|
|
_debug('But cond is always False.')
|
|
else:
|
|
_debug('Result before branch substitutions is:', res)
|
|
return _my_unpolarify(hyperexpand(res)), cond
|
|
|
|
# Try two G functions.
|
|
gs = _rewrite2(f, x)
|
|
if gs is not None:
|
|
for full_pb in [False, True]:
|
|
fac, po, g1, g2, cond = gs
|
|
_debug('Could rewrite as two G functions:', fac, po, g1, g2)
|
|
res = S.Zero
|
|
for C1, s1, f1 in g1:
|
|
for C2, s2, f2 in g2:
|
|
r = _rewrite_saxena(fac*C1*C2, po*x**(s1 + s2),
|
|
f1, f2, x, full_pb)
|
|
if r is None:
|
|
_debug('Non-rational exponents.')
|
|
return
|
|
C, f1_, f2_ = r
|
|
_debug('Saxena subst for yielded:', C, f1_, f2_)
|
|
cond = And(cond, _check_antecedents(f1_, f2_, x))
|
|
if cond == False:
|
|
break
|
|
res += C*_int0oo(f1_, f2_, x)
|
|
else:
|
|
continue
|
|
break
|
|
cond = _my_unpolarify(cond)
|
|
if cond == False:
|
|
_debugf('But cond is always False (full_pb=%s).', full_pb)
|
|
else:
|
|
_debugf('Result before branch substitutions is: %s', (res, ))
|
|
if only_double:
|
|
return res, cond
|
|
return _my_unpolarify(hyperexpand(res)), cond
|
|
|
|
|
|
def meijerint_inversion(f, x, t):
|
|
r"""
|
|
Compute the inverse laplace transform
|
|
$\int_{c+i\infty}^{c-i\infty} f(x) e^{tx}\, dx$,
|
|
for real c larger than the real part of all singularities of ``f``.
|
|
|
|
Note that ``t`` is always assumed real and positive.
|
|
|
|
Return None if the integral does not exist or could not be evaluated.
|
|
|
|
Examples
|
|
========
|
|
|
|
>>> from sympy.abc import x, t
|
|
>>> from sympy.integrals.meijerint import meijerint_inversion
|
|
>>> meijerint_inversion(1/x, x, t)
|
|
Heaviside(t)
|
|
"""
|
|
f_ = f
|
|
t_ = t
|
|
t = Dummy('t', polar=True) # We don't want sqrt(t**2) = abs(t) etc
|
|
f = f.subs(t_, t)
|
|
_debug('Laplace-inverting', f)
|
|
if not _is_analytic(f, x):
|
|
_debug('But expression is not analytic.')
|
|
return None
|
|
# Exponentials correspond to shifts; we filter them out and then
|
|
# shift the result later. If we are given an Add this will not
|
|
# work, but the calling code will take care of that.
|
|
shift = S.Zero
|
|
|
|
if f.is_Mul:
|
|
args = list(f.args)
|
|
elif isinstance(f, exp):
|
|
args = [f]
|
|
else:
|
|
args = None
|
|
|
|
if args:
|
|
newargs = []
|
|
exponentials = []
|
|
while args:
|
|
arg = args.pop()
|
|
if isinstance(arg, exp):
|
|
arg2 = expand(arg)
|
|
if arg2.is_Mul:
|
|
args += arg2.args
|
|
continue
|
|
try:
|
|
a, b = _get_coeff_exp(arg.args[0], x)
|
|
except _CoeffExpValueError:
|
|
b = 0
|
|
if b == 1:
|
|
exponentials.append(a)
|
|
else:
|
|
newargs.append(arg)
|
|
elif arg.is_Pow:
|
|
arg2 = expand(arg)
|
|
if arg2.is_Mul:
|
|
args += arg2.args
|
|
continue
|
|
if x not in arg.base.free_symbols:
|
|
try:
|
|
a, b = _get_coeff_exp(arg.exp, x)
|
|
except _CoeffExpValueError:
|
|
b = 0
|
|
if b == 1:
|
|
exponentials.append(a*log(arg.base))
|
|
newargs.append(arg)
|
|
else:
|
|
newargs.append(arg)
|
|
shift = Add(*exponentials)
|
|
f = Mul(*newargs)
|
|
|
|
if x not in f.free_symbols:
|
|
_debug('Expression consists of constant and exp shift:', f, shift)
|
|
cond = Eq(im(shift), 0)
|
|
if cond == False:
|
|
_debug('but shift is nonreal, cannot be a Laplace transform')
|
|
return None
|
|
res = f*DiracDelta(t + shift)
|
|
_debug('Result is a delta function, possibly conditional:', res, cond)
|
|
# cond is True or Eq
|
|
return Piecewise((res.subs(t, t_), cond))
|
|
|
|
gs = _rewrite1(f, x)
|
|
if gs is not None:
|
|
fac, po, g, cond = gs
|
|
_debug('Could rewrite as single G function:', fac, po, g)
|
|
res = S.Zero
|
|
for C, s, f in g:
|
|
C, f = _rewrite_inversion(fac*C, po*x**s, f, x)
|
|
res += C*_int_inversion(f, x, t)
|
|
cond = And(cond, _check_antecedents_inversion(f, x))
|
|
if cond == False:
|
|
break
|
|
cond = _my_unpolarify(cond)
|
|
if cond == False:
|
|
_debug('But cond is always False.')
|
|
else:
|
|
_debug('Result before branch substitution:', res)
|
|
from sympy.simplify import hyperexpand
|
|
res = _my_unpolarify(hyperexpand(res))
|
|
if not res.has(Heaviside):
|
|
res *= Heaviside(t)
|
|
res = res.subs(t, t + shift)
|
|
if not isinstance(cond, bool):
|
|
cond = cond.subs(t, t + shift)
|
|
from .transforms import InverseLaplaceTransform
|
|
return Piecewise((res.subs(t, t_), cond),
|
|
(InverseLaplaceTransform(f_.subs(t, t_), x, t_, None), True))
|