Rewriting expressions as a Meijer G-function for integration

[Ondřej Čertík]

Hi,

Our Meijer G-function integrator can find integrals (definite and
indefinite) of any function that can be expressed as a Meijer
G-function thanks to the formulas here:

http://functions.wolfram.com/HypergeometricFunctions/MeijerG/21/ShowAll.html

I.e. an integral of a Meijer G-function is also a Meijer G-function.
The definite integral has tons of conditions that SymPy checks, but
the formula for the indefinite integral (i.e. the antiderivative)
always holds.

Then one converts the final Meijer G-function into elementary
functions if possible, or leaves it as is if it is not possible. This
part is robust.

What is not robust is how to rewrite a given function into a Meijer
G-function. This is done by the `meijerint._rewrite1` function (btw,
we should expose it as a public function). For example:

In [1]: meijerint._rewrite1((cos(x)/x), x)
Out[1]: (1, 1/x, [(sqrt(pi), 0, meijerg(((), ()), ((0,), (1/2,)),
x**2/4))], True)

In [2]: meijerint._rewrite1((sin(x)/x), x)
Out[2]: (1, 1/x, [(sqrt(pi), 0, meijerg(((), ()), ((1/2,), (0,)),
x**2/4))], True)

In [3]: meijerint._rewrite1((cos(x)/x)**2, x)

In [4]: meijerint._rewrite1((sin(x)/x)**2, x)
Out[4]:
(1,
x**(-2),
[(sqrt(pi)/2, 0, meijerg(((0,), (1/2, 1/2, 1)), ((0, 1/2), ()), x**(-2)))],
True)

In [3] it didn't find the solution, yet a similar expression involving
sin(x) instead of cos(x) works in [4].

Let's figure out a systematic algorithm. For that, you start with the
elementary functions, that would be sin(x), cos(x) and "x" in the
above expression, look their G-function representation, and then use
the *, /, +, - and ** operations on the G-functions, that's it.

Now the problem is, that there doesn't seem to be a formula for a
product of two G functions, e.g. I didn't see one here:

http://functions.wolfram.com/HypergeometricFunctions/MeijerG/16/ShowAll.html

the formula http://functions.wolfram.com/HypergeometricFunctions/MeijerG/16/02/01/0001/
that you see there only seems to be using some even more generalized G
function of two variables? It doesn't seem to be useful here. Can
someone confirm that one cannot express a product of two G-functions
as a G-function?

So a solution is to simply have a robust method to rewrite any
expression as a hypergeometric function and then use the formula here
to convert the hypergeometric function to a G-function:

http://functions.wolfram.com/HypergeometricFunctions/HypergeometricPFQ/26/03/01/0001/

There are just a few functions that can be expressed as a G-function
but not as a hypergeometric function, some examples are: Bessel
functions Y, K (for integer order), Whittaker function W, Legendre
function Q_mu_nu and a few others. So for these functions we have to
figure out something else, probably something that we do now.

Also we can then use the integration formulas here for hypergeometric
functions, so we don't even have to go via G-functions:

http://functions.wolfram.com/HypergeometricFunctions/HypergeometricPFQ/21/ShowAll.html

It seems the conditions on definite integration are a lot simpler as well.


So here is the algorithm for hypergeometric functions, I'll show it on
the (sin(x)/x)**2 example above:

1) sin(x) = x * 0F1(3/2; -x^2/4)

2) sin(x) / x = 0F1(3/2; -x^2/4)

3) (sin(x)/x)**2 = 0F1(3/2; -x^2/4) * 0F1(3/2; -x^2/4) = 2F3(3/2,1;
3/2,3/2,2; -x^2)

Where we used the formula for a product of two 0F1 functions:

http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric0F1/16/ShowAll.html

4) Finally, rewrite 2F3(3/2,1; 3/2,3/2,2; -x^2) as a G-function, or
integrate directly.


A general formula for multiplication of two hypergeometric series is here:

http://functions.wolfram.com/HypergeometricFunctions/HypergeometricPFQ/16/ShowAll.html

But I can see now that this only expresses the result as a Taylor
series. So maybe I just got lucky that the multiplication of two 0F1
functions exists (in terms of 2F3), but there is no such formula for a
general PFQ function. Can someone confirm this?

So then the other idea is that one can identify a hypergeometric
function from the series expansion. Essentially, we expand into a
series, calculate the ratio t_{k+1}/t_k of two successive terms, and
if it is a rational function of "k", then it is a hypergeometric
function and you read the coefficients directly. Let's do the same
example again:

1) sin^2(x) = Sum(((-1)**(k - 1)* 2**(2* k - 1)*
x**(2*k))/factorial(2*k), (k, 1, oo))

2) (sin(x)/x)**2 = sin^2(x)/x^2 = Sum(((-1)**(k - 1)* 2**(2* k - 1)*
x**(2*k-2))/factorial(2*k), (k, 1, oo))
= Sum(((-1)**k * 2**(2*k+1) * x**(2*k))/factorial(2*k+2), (k, 0, oo))
= Sum((2**(2*k+1) * (-x**2)**k))/factorial(2*k+2), (k, 0, oo))

So the series is of the form sum_k t_k*(-x^2)^k where

t_k = 2**(2*k+1)/factorial(2*k+2)

3) Calculate t_{k+1} / t_k:

In [78]: t_k = 2**(2*k+1)/factorial(2*k+2)

In [79]: t_k.subs(k, k+1) / t_k
Out[79]: 2**(-2*k - 1)*2**(2*k + 3)*factorial(2*k + 2)/factorial(2*k + 4)

In [80]: _.simplify()
Out[80]: 2/((k + 2)*(2*k + 3))

We write the last formula as:

2/((k + 2)*(2*k + 3)) = (k+1) / ((k+3/2)*(k+2)*(k+1))

And we read off the hypergeometric function as 1F2(1; 3/2, 2; -x^2).

One can verify that we got the same answer as before:

In [64]: hyperexpand(hyper([1], [S(3)/2, 2], -x**2))
Out[64]:
cos(2⋅x) 1
- ──────── + ────
2 2
2⋅x 2⋅x

In [65]: hyperexpand(hyper([S(3)/2, 1], [S(3)/2, S(3)/2, 2], -x**2))
Out[65]:
cos(2⋅x) 1
- ──────── + ────
2 2
2⋅x 2⋅x


And both are actually equal to (sin(x)/x)^2:

In [75]: (_ - (sin(x)/x)**2).simplify()
Out[75]: 0


One can see, that one should use the exact series expansion, what is
the status of that, I think we had a GSoC on it?


It might be, that all the operations that I am doing on the series
have an equivalent operation on the hypergeometric series.

In conclusion, the logic is simple: either the product or power of
hypergeometric functions exists as a hypergeometric function for the
given set of coefficients (like above) and then there must be an
automated way to determine the final hypergeometric function (in the
worst case by calculating the t_{k+1}/tk ratio like I did above), or
the answer can't be expressed as a hypergeometric function (i.e. when
the t_{k+1}/tk ratio is not a rational function of "k"), and then the
answer is that there is no hypergeometric function that represents the
result, and so we can't integrate it using this method and that's
fine.

Let me know what you think.


This was motivated by https://github.com/sympy/sympy/issues/10464.

Ondrej

[Ondřej Čertík]

Btw, this procedure is explained at the page 36 of the freely
available A=B book, together with many examples:

https://www.math.upenn.edu/~wilf/Downld.html

As they show, one can definitely determine if a given series is
hypergeometric (and find the function) or not. So the only hard part
is to find the series expansion of the final expression in closed form
(i.e. as an infinite sum, but have a closed form expression for each
coefficient) and then apply this algorithm.

Looking at the chapters 4 and 5, I think we can actually make use of
the general multiplication formula here:

http://functions.wolfram.com/HypergeometricFunctions/HypergeometricPFQ/16/ShowAll.html

Since it takes two general hypergeometric functions, multiplies them
and writes the result as an infinite series: Sum(c_k * z**k, (k, 0,
oo)), where the coefficients c_k are given using a hypergeometric
function. Then we apply the algorithm from page 36 to determine
whether or not this sum can be written as a hypergeometric function,
and if so, which one (i.e. by calculating and simplifying
c_{k+1}/c_k). So this is what we need to implement. And the chapters 4
and 5 also treat similar sums, where the c_k coefficients are given
using a hypergeometric function.

Ondrej

[Ondřej Čertík]

Here is a script that implements it:

https://gist.github.com/certik/a47e2040d9a44812e2f2

I tested a few examples (see the commented out lines in the script)
and it works as expected.

One starts with two hypergeometric sequences, it calculates the
coefficients in the infinite series of the product, then calculates
the ratio. E.g. for the example above the ratio is:

-gamma(k + 3/2)/((k + 2)*gamma(k + 5/2))

It would be nice if sympy could simplify it to:

-1/((k + 2)*(k+3/2)) = -(k+1)/((k + 2)*(k+3/2)*(k+1))

But in the meantime I simply read the hypergeometric function manually
in the script as hyper([1], [2, S(3)/2], -z) in this case. Then the
script compares the series expansions to show that it works.

So I think this is a viable approach.

Ondrej



P.S. One thing that I don't understand is the cos^2(x) case. The
general formula is:

0F1^2(b; z) = 2F3(b, b-1/2; b, b, 2b-1; 4z) = 1F2(b-1/2; b, 2b-1; 4z)

E.g. for sin(x)/x we have b=3/2 and get 0F1^2(3/2; -x^2/4) = 1F2(1;
3/2, 2; -x^2), consistent with above. But for cos(x) we get b=1/2 and:

cos^2(x) = 0F1^2(1/2; -x^2/4) = 1F2(0; 1/2, 0; -x^2)

Which is not well defined. However, I have experimentally found out, that:

cos^2(x) = 0F1(1/2; -x^2) + x^2*1F2(1; 3/2, 2; -x^2)

>>> hyperexpand(hyper([], [S(1)/2], -x**2) + x**2*hyper([1], [S(3)/2, 2], -x**2)).trigsimp()
cos(x)**2

The script for this case returns a ratio of:

-gamma(k + 1/2)/((k + 1)*gamma(k + 3/2)) = -1/((k+1)*(k+1/2))

Which would suggest 0F1(1/2; -x^2) as the solution, but that is incorrect:

>>> hyperexpand(hyper([], [S(1)/2], -x**2))
cos(2*x)

since cos(2*x) = cos^2(x) - sin^2(x) and one needs to add
sin^2(x)=x^2*1F2(1; 3/2, 2; -x^2) to it, as shown above. Since this
case is essentially a linear combination of hypergeometric functions,
maybe something in the algorithm fails. I would expect the ratio to
return something that is not a rational function, since the result is
not just one hypergeometric function. One would have to look into this
more.

[Aaron Meurer]

I'm CCing Tom Bachmann, since he wrote the original code.

+1 for exposing a function that can rewrite expressions to hyper or meijerg as a public function (perhaps just expr.rewrite(hyper) and expr.rewrite(meijerg)).

Aaron Meurer

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To post to this group, send email to sy...@googlegroups.com.
Visit this group at https://groups.google.com/group/sympy.
To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CADDwiVBm6KhcTOkf8Smo86Ucro5CAx6LOLeCRpzYZ5wNu8paqg%40mail.gmail.com.

For more options, visit https://groups.google.com/d/optout.

[Kalevi Suominen]

To get 0F1^2(1/2; x) from the formula

0F1^2(b; z) = 1F2(b-1/2; b, 2b-1; 4z)

it seems necessary to consider the limit  b -> 1/2. The coefficients of the 1F2 series,
after the constant term 1,

rf(b-1/2, n)/rf(b, n)rf(2b-1, n)n! = (b-1/2)rf(b+1/2, n-1)/(2b-1)rf(2b, n-1)rf(b, n)n!

become

rf(1, n-1)/2rf(1, n-1)rf(1/2, n)n! = 1/2rf(1/2, n)n!,

i.e., one half of the coefficients of 0F1(1/2; z).
Hence we get

0F1^2(1/2; z) = 1 + (0F1(1/2; 4z) - 1)/2 = (1 + 0F1(1/2; 4z))/2

which is not a hypergeometric series (as a series of 4z).
The ratio of its consecutive coefficients is the same as
in  0F1 except at the first step where it is only one half
of the expected value.

Kalevi

[Ondřej Čertík]

Hi Kalevi,

You nailed it! I checked your calculations and they are correct. Your
last equation is nothing else than:

cos(x)^2 = 0F1^2(1/2; -x^2/4) = (1+0F1(1/2; -x^2))/2 = (1+cos(2*x))/2

Which is obviously true. So for integration purposes, the final answer
(1+0F1(1/2; -x^2))/2 is fine, as it is a linear combination of
hypergeometric functions, which can be integrated easily. You are
right, that cos^2(x) is not a (single) hypergeometric series. Which is
fine, there is problem.

Thanks for doing this, it was bothering me that things didn't work for
b=1/2. They work just fine.

Ondrej

[Ondřej Čertík]

On Mon, Feb 1, 2016 at 1:26 PM, Ondřej Čertík <ondrej...@gmail.com> wrote:
[...]

> right, that cos^2(x) is not a (single) hypergeometric series. Which is
> fine, there is problem.

-> there is no problem.

Ondrej

[brandon willard]

I've been thinking about this same topic a lot recently (partially due to a question about a G-function form of tanh), and it seems like the more generalized G-function you mentioned, Ondrej, is probably necessary at some point.  There doesn't seem to be a whole lot of literature on these bivariate G-functions, but, if you extend the scope to H-functions and bivariate hypergeometric functions (e.g. Horn, Appell), there are at least enough useful identities to consider implementing.

Here's one interesting identity involving that generalized G-function and the Appell: http://functions.wolfram.com/07.34.16.0003.01.

Also, there are some explicit series expansions for H-functions might help: http://arxiv.org/abs/math/9803163.

[Ondřej Čertík]

Hi,

I just want to update this thread, that we have fixed this problem, in
my opinion, thanks to our GSoC student Subham Tibra, and his mentor
Kalevi Suominen. I am Subham's mentor too, but Kalevi has done a much
better job mentoring. ;)

In short: hypergeometric as well as the MeijerG functions are
solutions to an ODE together with some (symbolic) initial conditions
at a point like x=0, or x=1 (or any other point). Holonomic functions
are a generalization of this ODE to allow polynomials as coefficients
of the ODE. Then suddenly they are closed to almost all operations
(including multiplication, integration, differentiation) and there are
algorithms that can robustly and quickly compute those operations.
They are now implemented in SymPy.

I am using SymPy version 8d7b522e58aae883b4592e4fae3babf82d1e4db2.
Let's first show that what the meijerint._rewrite1() cannot do, can be
done easily with holonomic functions:

In [1]: from sympy.holonomic import from_sympy

In [2]: meijerint._rewrite1((cos(x)/x), x)
Out[2]: (1, 1/x, [(sqrt(pi), 0, meijerg(((), ()), ((0,), (1/2,)),
x**2/4))], True)

In [3]: meijerint._rewrite1((sin(x)/x), x)
Out[3]: (1, 1/x, [(sqrt(pi), 0, meijerg(((), ()), ((1/2,), (0,)),
x**2/4))], True)

In [4]: meijerint._rewrite1((cos(x)/x)**2, x)

In [5]: meijerint._rewrite1((sin(x)/x)**2, x)
Out[5]:

(1,
x**(-2),
[(sqrt(pi)/2, 0, meijerg(((0,), (1/2, 1/2, 1)), ((0, 1/2), ()), x**(-2)))],
True)

In [6]: from_sympy((cos(x)/x))
Out[6]: HolonomicFunction((x) + (2)Dx + (x)Dx**2, x), f(1) = cos(1),
f'(1) = -sin(1) - cos(1)

In [7]: from_sympy((sin(x)/x))
Out[7]: HolonomicFunction((x) + (2)Dx + (x)Dx**2, x), f(0) = 1, f'(0) = 0

In [8]: from_sympy((cos(x)/x)**2)
Out[8]: HolonomicFunction((8*x) + (4*x**2 + 6)Dx + (6*x)Dx**2 +
(x**2)Dx**3, x), f(1) = cos(1)**2, f'(1) = -2*sin(1)*cos(1) -
2*cos(1)**2, f''(1) = 4*cos(1)**2 + 2*sin(1)**2 + 8*sin(1)*cos(1)

In [9]: from_sympy((sin(x)/x)**2)
Out[9]: HolonomicFunction((8*x) + (4*x**2 + 6)Dx + (6*x)Dx**2 +
(x**2)Dx**3, x), f(0) = 1, f'(0) = 0, f''(0) = -2/3



Here is an example how to use the holonomic functions module to
compute integrals:

https://github.com/sympy/sympy/issues/8944#issuecomment-229478358

it's around 10x faster than the SymPy's integrate() routine.

The tough part is what to do about definite integrals where the
antiderivative (a holonomic function) can be converted to elementary
functions, like here:

https://github.com/sympy/sympy/issues/11319

That's where the MeijerG approach gives better results. Also another
advantage of the MeijerG approach is that it gives convergence
conditions --- though perhaps there is a way to implement it in the
holonomic module (https://github.com/sympy/sympy/issues/11322).

Ondrej

P.S. Thanks Brandon for your email. I think the above is the solution.

> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at https://groups.google.com/group/sympy.
> To view this discussion on the web visit

> https://groups.google.com/d/msgid/sympy/3e750003-cf00-4e6e-8e1c-c032fc7dc037%40googlegroups.com.

[brandon willard]

Brilliant! I have come across the holonomic/[Meijer] G-function connection before, but never really explored it. Do G/H-functions represent the entire space of holonomic functions? How about the other way around?
I also recall this thesis covering both in connection to the solution of integrals.

​

You received this message because you are subscribed to a topic in the Google Groups "sympy" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/sympy/S0MSlsAhL74/unsubscribe.
To unsubscribe from this group and all its topics, send an email to sympy+un...@googlegroups.com.

To post to this group, send email to sy...@googlegroups.com.
Visit this group at https://groups.google.com/group/sympy.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CADDwiVBvFqR6FF3cdn%3Dh9r57xaJk70JM1yQ_M2E91DZ5Rh%2BK0g%40mail.gmail.com.

[Ondřej Čertík]

Hi Brandon,

Thanks for the link. No, G-functions are a subset of holonomic
functions. In some sense, the G-functions are not complete, since they
are not closed under multiplication. But holonomic functions are.
Thanks for the link to the thesis, Subham, you should have a look at
it, it might contains some useful algorithms.

You can follow our progress on the holonomic package by following our issues:

https://github.com/sympy/sympy/labels/Holonomic%20Functions

One problem that I didn't realize before is this:
https://github.com/sympy/sympy/issues/11329, that from the holonomic
differential equation, it's easy to identify, if the equation is an
equation for a G-function, but the problem is that G-functions are
various solutions of the corresponding ODE, so one has to identify the
right linear combination, based on the initial condition.

A similar problem is when converting to hypergeometric functions. But
there we use the fact that we can obtain a closed recursive formula
for the coefficients of a Taylor series of the given holonomic
function (from that one can identify the hypergeometric function).

Ondrej

On Thu, Jun 30, 2016 at 2:17 PM, brandon willard

> https://groups.google.com/d/msgid/sympy/CAEOXDQ2HM9a%3DUUGnk6%2BAZZOnpBpMhyz2Ki-BCEn7oa7qQJKa6A%40mail.gmail.com.

[Raymond Rogers]

"Can someone confirm that one cannot express a product of two G-functions as a G-function? "

Roundabout but:

If one has the G-function expression for the Mellin transform of the convolution of two functions then the inverse Mellin operator applied to it would yield the G-function for the product?

I am not an expert in sympy or G-functions; just working my way through

"The Double Mellin-Barnes Type Integrals and Their Applications to Convolution Theory"

and sometimes wandering around in my mind.

RayR