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 On Wed, Jan 27, 2016 at 10:14 PM, Ondřej Čertík <ondrej.cer...@gmail.com> wrote: > On Wed, Jan 27, 2016 at 4:30 PM, Ondřej Čertík <ondrej.cer...@gmail.com> > wrote: > > On Wed, Jan 27, 2016 at 3:34 PM, Ondřej Čertík <ondrej.cer...@gmail.com> > wrote: > >> 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). > > > > 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. > > 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. > > -- > 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+unsubscr...@googlegroups.com. > To post to this group, send email to sympy@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. > -- 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+unsubscr...@googlegroups.com. To post to this group, send email to sympy@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/CAKgW%3D6%2Ba6T2E5W1OSTeudkSfCEOAfhJCkvYG_hGFK-6A6eNR4w%40mail.gmail.com. For more options, visit https://groups.google.com/d/optout.
