How Mathematica and Maple Get Meijer's G-function into Problem 9

Folkmar Bornemann, Feb 8, 2005

in memoriam John Boersma (1937-2004)

Asking Wolfram Research, Inc.

Stan's inquiry to the WRI team as to how Mathematica manages to get Meijer's G-function into play, was answered by Daniel Lichtblau on July 27, 2003,  in the following rather cryptic way:

I'll be discussing aspects of MeijerG issues at ACA next week. It is basically a lookup that converts various functions to MeijerG, then figures out the integral of a product of 2 MeijerG's via Slater convolution.

This was not detailed enough to help our understanding of how this works.

Mathematica Reference Guide

The Notes on Internal Implementation reveal the following knowledge:

Many other definite integrals are done using Marichev-Adamchik Mellin transform methods. The results are often initially expressed in terms of Meijer G functions, which are converted into hypergeometric functions using Slater's Theorem and then simplified.

The problem is similar as above: these keywords are not sufficient to get reasonable pointers to some literature. But there are some names mentioned which could be traced.

The Relevant Literature

An extensive search with Google Scholar (great stuff indeed) has led me to the following relevant items:

[Mari83] Oleg Marichev, Handbook of Integral Transforms of Higher Transcendental Functions: Theory and Algorithmic Tables, Ellis Horwood Ltd., Chichester, 1983.

[Adam96] Victor Adamchik, Definite Integration in Mathematica,  Mathematica in Education and Research, no. 3, 5 (1996), 16–22.

According to [Adam96] the ideas of [Mari83] are implemented in Mathematica. Marichev's involvement with Maple suggests that the same is true for this CAS, too.

The General Idea

Marichev's method is based on the observation that any definite integral can be written in the form of Mellin's convolution of two functions, that is

f (z)=(∫_0^∞)(Subscript[f, 1])(x)(Subscript[f, 2])((z/x))d(x/x).

Now, this is a convolution since the Mellin transforms of these functions, that is,

\(\(f\^*\) \((s)\) = ∫\_0\%∞ 
            f \((x)\) \(x\^\(s - 
                  1\)\) dx, \ \ \ \ \ \ \ \ \ \ \ \(f\_j\%*\) \((s)\) = ∫\_0\%∞\( 
                f\_j\) \((x)\) \(x\^\(s - 
                  1\)\) dx, \ \ \ \ \ j = 1, 2, \)

simply multiply

\(\(f\^*\) \((s)\) = \ \(f\_1\%*\) \((s)\)·\ \(f\_2\%*\) \(\((s)\)\(.\)\)\)

Of course, generally this observation will be of not much use. However, if the functions are of hypergeometric type, which is true for many elementary functions and the majority of special functions, the inverse Mellin transform of f,

\(\(\(f \((z)\) = \(1\/\(2  π\ i\)\) \(∳\(f\^*\) \((s)\) \(z\^\(\(-s\)\(\ \)\)\) ds\)\)\(,\)\)\)

turns out to be a Mellin-Barnes integral. (One has to exercise some care in choosing a path of integration.) Depending on the involved coefficients, this may represent Fox's H-function, or in simpler cases, Meijer's G-function. Finally, Lucy Slater's theorem of 1966 (see [Mari83, p. 57]) states precise conditions on those coefficients such that Meijer's G-function reduces to a hypergeometric series, and therefore, in many cases can be simplified to more familiar special functions. In fact, this situation of a reduction to a hypergeometric series turns out to be the generic case.

In the non-reducible, non-generic case, one can still evaluate the Mellin-Barnes integral by means of the residue theorem yielding an infinite series involving polygamma functions. This way Maple 9.5 completely avoids delivering Meijer's G-functions.

Let's Do It: Problem 9

We want to use this information to understand what Mathematica is doing when performing the following evaluation:

J = 
      Simplify[(∫_0^2)(((2 - x))^α) sin((α/x)) \[DifferentialD]x, {α > 0}]

With z = 1 and

(Subscript[f, 1])(x_) = 
      x (((2 - x))^α) If[0 ≤ x ≤ 2, 1, 
          0];            (Subscript[f, 2])(x_) = 
      sin(α x);

the integral is obviously of the required convolution type. The Mellin transforms of (Subscript[f, 1]) and (Subscript[f, 2]) are given by

\(TraditionalForm\`\(f\_1\%*\)(s_) = 
      Simplify[∫\_0\%∞\(\( f\_1\)(
              x)\)\ \(x\^\(s - 1\)\) \[DifferentialD]x, {α > 0, 
          Re(s) > \(-1\)}]\)

and—by the reflection formula of the gamma function—

\(TraditionalForm\`\(f\_2\%*\)(s_) = 
      Simplify[∫\_0\%∞\(\( f\_2\)(
                x)\)\ \(x\^\(s - 
                    1\)\) \[DifferentialD]x /. \((sin(π\ t_) -> \ π\/\(\(Γ(t)\) \(Γ(1 - t)\)\))\), {α > 
            0, 0 < Re(s) < 1}]\)

Both transforms, which could also have been obtained by a short calculation by hand, are of the required hypergeometric form.

Thus, the inverse Mellin transform yields the Mellin-Barnes integral

J=f (1)=(1/(2  π i))∳((π (2^(α + 1))Γ(α + 1)Γ(s + 1) Γ(s))/(Γ(s + α + 2)Γ(s/2) Γ(
                        1 - s/2)))((((α/2)))^(-s))ds .

This is not yet recognizable as a Meijer's G-function, since in the definition of these functions the variable s  can only enter the gamma function terms with a coefficient ±1. Therefore we apply the transform s ->  2  s and obtain—by using the duplication formula for the gamma function—the transformed integrand as:

\(TraditionalForm\`\(2 \(\( f\_1\%*\)(s)\) \(\(f\_2\%*\)(
              s)\) /. \((s -> 
              2\ s)\)\) /. \((Γ(
              t : 2  γ_ + δ_:  0) -> \(1\/\@π\) 
              2\^\(t - 1\)\ \(Γ(
                t\/2)\)\ \(Γ(\(t + 1\)\/2)\))\) // 
      Simplify\)

Summarizing we get

J=f(1)=(Sqrt[π])Γ(α + 1)·((1/(2  π i))∳((((Γ(
                                      s + (1/2)))^2) Γ(
                                  s + 1) )/(\(\(Γ(
                                      s + \(α + 2\)\/2)\)\ \(Γ(s + \(α + 3\)\/2)\)\)Γ(1 - s) ))(((((α^2)/16)))^(-s))ds ),

where the large parenthesis is nothing but the definition of Meijer's G-function,

J=(Sqrt[π])Γ(α + 1)·(Subscript[G, 2, 4]^(3, 0))(((α^2)/16)|({{((α + 2)/2),((α + 3)/2)}, {(1/2),(1/2),1,0}})).

(It does not reduce to a hypergeometric series because of the double appearance of the coefficient (1/2).)

Thus, the approach is very systematic while involving a minimum of calculation.