Computing Integrals using Meijer G-Functions¶
This text aims do describe in some detail the steps (and subtleties) involved in using Meijer G-functions for computing definite and indefinite integrals. We shall ignore proofs completely.
Overview¶
The algorithm to compute
Rewrite the integrand using Meijer G-functions (one or sometimes two).
Apply an integration theorem, to get the answer (usually expressed as another G-function).
Expand the result in named special functions.
Step (3) is implemented in the function hyperexpand (q.v.). Steps (1) and (2) are described below. Moreover, G-functions are usually branched. Thus our treatment of branched functions is described first.
Some other integrals (e.g.
Polar Numbers and Branched Functions¶
Both Meijer G-Functions and Hypergeometric functions are typically branched
(possible branchpoints being
To begin, we consider the set
We also sometimes formally attach a point “zero” (
Using these maps many operations can be extended from
As one peculiarity it should be mentioned that addition of polar numbers is not
usually defined. However, formal sums of polar numbers can be used to express
branching behaviour. For example, consider the functions
(We are omitting
Representing Branched Functions on the Argand Plane¶
Suppose
Table Lookups and Inverse Mellin Transforms¶
Suppose we are given an integrand _split_mul(f, x)
. Then we assemble a tuple of functions that occur in
_mytype(f, x)
. Next we index a lookup table
(created using _create_lookup_table()
) with this tuple. This (hopefully)
yields a list of Meijer G-function formulae involving these functions, we then
pattern-match all of them. If one fits, we were successful, otherwise not and we
have to try something else.
Suppose now we want to rewrite as a product of two G-functions. To do this,
we (try to) find all inequivalent ways of splitting _mul_as_two_parts(f)
.
Finally, we can try a recursive Mellin transform technique. Since the Meijer
G-function is defined essentially as a certain inverse mellin transform,
if we want to write a function
One twist is that some functions don’t have mellin transforms, even though they
can be written as G-functions. This is true for example for
The function _rewrite_single
does the table lookup and recursive mellin
transform. The functions _rewrite1
and _rewrite2
respectively use
above-mentioned helpers and _rewrite_single
to rewrite their argument as
respectively one or two G-functions.
Applying the Integral Theorems¶
If the integrand has been recast into G-functions, evaluating the integral is
relatively easy. We first do some substitutions to reduce e.g. the exponent
of the argument of the G-function to unity (see _rewrite_saxena_1
and
_rewrite_saxena
, respectively, for one or two G-functions). Next we go through
a list of conditions under which the integral theorem applies. It can fail for
basically two reasons: either the integral does not exist, or the manipulations
in deriving the theorem may not be allowed (for more details, see this [BlogPost]).
Sometimes this can be remedied by reducing the argument of the G-functions
involved. For example it is clear that the G-function representing meijerg.get_period()
can be used to discover this, and the function
principal_branch(z, period)
in functions/elementary/complexes.py
can
be used to exploit the information. This is done transparently by the
integration code.
The G-Function Integration Theorems¶
This section intends to display in detail the definite integration theorems used in the code. The following two formulae go back to Meijer (In fact he proved more general formulae; indeed in the literature formulae are usually staded in more general form. However it is very easy to deduce the general formulae from the ones we give here. It seemed best to keep the theorems as simple as possible, since they are very complicated anyway.):
The more interesting question is under what conditions these formulae are valid. Below we detail the conditions implemented in SymPy. They are an amalgamation of conditions found in [Prudnikov1990] and [Luke1969]; please let us know if you find any errors.
Conditions of Convergence for Integral (1)¶
We can without loss of generality assume
The convergence conditions will be detailed in several “cases”, numbered one
to five. For later use it will be helpful to separate conditions “at infinity”
from conditions “at zero”. By conditions “at infinity” we mean conditions that
only depend on the behaviour of the integrand for large, positive values
of
In order for the integral theorem to be valid, conditions
These are the conditions “at infinity”:
where
And these are the conditions “at zero”:
Conditions of Convergence for Integral (2)¶
We introduce the following notation:
The following conditions will be helpful:
Note that
With this notation established, the implemented convergence conditions can be enumerated as follows:
The Inverse Laplace Transform of a G-function¶
The inverse laplace transform of a Meijer G-function can be expressed as
another G-function. This is a fairly versatile method for computing this
transform. However, I could not find the details in the literature, so I work
them out here. In [Luke1969], section 5.6.3, there is a formula for the inverse
Laplace transform of a G-function of argument
We are asked to compute
for positive real
When does this integral converge?
How can we compute the integral?
When is our computation valid?
How to compute the integral¶
We shall work formally for now. Denote by
Thus
We interchange the order of integration to get
The inner integral is easily seen to be
which is a so-called Fox H function (of argument
When this computation is valid¶
There are a number of obstacles in this computation. Interchange of integrals is only valid if all integrals involved are absolutely convergent. In particular the inner integral has to converge. Also, for our identification of the final integral as a Fox H / Meijer G-function to be correct, the poles of the newly obtained gamma function must be separated properly.
It is easy to check that the inner integral converges absolutely for
It remains to observe that the Meijer G-function is an analytic, unbranched
function of its parameters, and of the coefficient
When the integral exists¶
Using [Luke1969], section 5.10, for any given meijer G-function we can find a
dominant term of the form
We must thus investigate
(This principal value integral is the exact statement used in the Laplace
inversion theorem.) We write
Set
or . In this case the integral converges if . , , . In this case the integral converges if . , , , , and at least one of . Here the same condition as in (1) applies.
Implemented G-Function Formulae¶
An important part of the algorithm is a table expressing various functions as Meijer G-functions. This is essentially a table of Mellin Transforms in disguise. The following automatically generated table shows the formulae currently implemented in SymPy. An entry “generated” means that the corresponding G-function has a variable number of parameters. This table is intended to shrink in future, when the algorithm’s capabilities of deriving new formulae improve. Of course it has to grow whenever a new class of special functions is to be dealt with.
Elementary functions:
Functions involving
Functions involving
Functions involving
Functions involving
Functions involving
Functions involving
Functions involving
Functions involving
Functions involving
Functions involving
Functions involving
Functions involving
Functions involving
Functions involving
Functions involving
Functions involving
Functions involving
Functions involving
Functions involving
Functions involving
Functions involving
Functions involving
Functions involving
Functions involving
Functions involving
Functions involving
Functions involving
Internal API Reference¶
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
- exception sympy.integrals.meijerint._CoeffExpValueError[source]¶
Exception raised by _get_coeff_exp, for internal use only.
- sympy.integrals.meijerint._check_antecedents(g1, g2, x)[source]¶
Return a condition under which the integral theorem applies.
- sympy.integrals.meijerint._check_antecedents_1(g, x, helper=False)[source]¶
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
.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 exists.
- sympy.integrals.meijerint._check_antecedents_inversion(g, x)[source]¶
Check antecedents for the laplace inversion integral.
- sympy.integrals.meijerint._condsimp(cond, first=True)[source]¶
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
- sympy.integrals.meijerint._create_lookup_table(table)[source]¶
Add formulae for the function -> meijerg lookup table.
- sympy.integrals.meijerint._dummy(name, token, expr, **kwargs)[source]¶
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.
- sympy.integrals.meijerint._dummy_(name, token, **kwargs)[source]¶
Return a dummy associated to name and token. Same effect as declaring it globally.
- sympy.integrals.meijerint._exponents(expr, x)[source]¶
Find the exponents of
x
(not including zero) inexpr
.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}
- sympy.integrals.meijerint._find_splitting_points(expr, x)[source]¶
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}
- sympy.integrals.meijerint._flip_g(g)[source]¶
Turn the G function into one of inverse argument (i.e. G(1/x) -> G’(x))
- sympy.integrals.meijerint._functions(expr, x)[source]¶
Find the types of functions in expr, to estimate the complexity.
- sympy.integrals.meijerint._get_coeff_exp(expr, x)[source]¶
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)
- sympy.integrals.meijerint._guess_expansion(f, x)[source]¶
Try to guess sensible rewritings for integrand f(x).
- sympy.integrals.meijerint._inflate_fox_h(g, a)[source]¶
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.
- sympy.integrals.meijerint._inflate_g(g, n)[source]¶
Return C, h such that h is a G function of argument z**n and g = C*h.
- sympy.integrals.meijerint._int0oo(g1, g2, x)[source]¶
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(((0, 1/2), ()), ((m/2,), (-m/2,)), s**(-2))/s**2
- sympy.integrals.meijerint._int0oo_1(g, x)[source]¶
Evaluate
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))
- sympy.integrals.meijerint._int_inversion(g, x, t)[source]¶
Compute the laplace inversion integral, assuming the formula applies.
- sympy.integrals.meijerint._is_analytic(f, x)[source]¶
Check if f(x), when expressed using G functions on the positive reals, will in fact agree with the G functions almost everywhere
- sympy.integrals.meijerint._meijerint_definite_2(f, x)[source]¶
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.
- sympy.integrals.meijerint._meijerint_definite_3(f, x)[source]¶
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.
- sympy.integrals.meijerint._meijerint_definite_4(f, x, only_double=False)[source]¶
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.
- sympy.integrals.meijerint._meijerint_indefinite_1(f, x)[source]¶
Helper that does not attempt any substitution.
- sympy.integrals.meijerint._mul_args(f)[source]¶
Return a list
L
such thatMul(*L) == f
.If
f
is not aMul
orPow
,L=[f]
. Iff=g**n
for an integern
,L=[g]*n
. Iff
is aMul
,L
comes from applying_mul_args
to all factors off
.
- sympy.integrals.meijerint._mul_as_two_parts(f)[source]¶
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))]
- sympy.integrals.meijerint._my_principal_branch(expr, period, full_pb=False)[source]¶
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.
- sympy.integrals.meijerint._mytype( ) tuple[type[Basic], ...] [source]¶
Create a hashable entity describing the type of f.
- sympy.integrals.meijerint._rewrite1(f, x, recursive=True)[source]¶
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 ofx
. and po = x**s. Here g is a result from _rewrite_single. Return None on failure.
- sympy.integrals.meijerint._rewrite2(f, x)[source]¶
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.
- sympy.integrals.meijerint._rewrite_saxena(fac, po, g1, g2, x, full_pb=False)[source]¶
Rewrite the integral
fac*po*g1*g2
from 0 to oo in terms of G functions with argumentc*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)
- sympy.integrals.meijerint._rewrite_saxena_1(fac, po, g, x)[source]¶
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.
- sympy.integrals.meijerint._rewrite_single(f, x, recursive=True)[source]¶
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.
- sympy.integrals.meijerint._split_mul(f, x)[source]¶
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))
- sympy.integrals.meijerint.meijerint_definite(f, x, a, b)[source]¶
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.
- sympy.integrals.meijerint.meijerint_indefinite(f, x)[source]¶
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)
- sympy.integrals.meijerint.meijerint_inversion(f, x, t)[source]¶
Compute the inverse laplace transform
, for real c larger than the real part of all singularities off
.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)