Hi Kalevi, On Mon, Feb 1, 2016 at 2:11 AM, Kalevi Suominen <[email protected]> wrote: > > > 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 <[email protected]> >> wrote: >> > On Wed, Jan 27, 2016 at 3:34 PM, Ondřej Čertík <[email protected]> >> > 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.
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 -- 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 [email protected]. To post to this group, send email to [email protected]. 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/CADDwiVDYM2oQEaWzVp7%3D6tb3cwquW4d31rifq4jpWXNYrqcQXA%40mail.gmail.com. For more options, visit https://groups.google.com/d/optout.
