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. Ondrej > > 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 -- 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/CADDwiVBLjAHCsVDB4_x%2B891neaBp%2BUZbboj5A1TVKf2pB%3DeLyA%40mail.gmail.com. For more options, visit https://groups.google.com/d/optout.
