On Thursday, January 28, 2016 at 5:14:41 AM UTC+2, Ondřej Čertík wrote: > > On Wed, Jan 27, 2016 at 4:30 PM, Ondřej Čertík <ondrej...@gmail.com > <javascript:>> wrote: > > On Wed, Jan 27, 2016 at 3:34 PM, Ondřej Čertík <ondrej...@gmail.com > <javascript:>> 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. >
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 -- 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/2e846e81-54b3-4dc1-8644-8ee9edc10446%40googlegroups.com. For more options, visit https://groups.google.com/d/optout.
