MathGroup Archive 2004

[Date Index] [Thread Index] [Author Index]

Search the Archive

Hypergeometric and MeijerG

  • To: mathgroup at smc.vnet.net
  • Subject: [mg49275] Hypergeometric and MeijerG
  • From: ab_def at prontomail.com (Maxim)
  • Date: Sat, 10 Jul 2004 02:48:57 -0400 (EDT)
  • Sender: owner-wri-mathgroup at wolfram.com

There are some issues with numerical evaluation of the hypergeometric
and MeijerG functions. Consider:

In[1]:=
Hypergeometric1F1[-n, -3*n, 2] /. n -> 3 + 1`20*^-20
Hypergeometric1F1[-n, -3*n, 2] /. n -> 3 + 1*^-20 // N[#, 20]&

Out[1]=
1.8492188014784114060740699412`19.46436426588571

Out[2]=
1.84920634920634920634920634920634920635`20.

The results are different, but both have only about six correct
decimal digits, which can be seen by setting the precision higher.
This kind of errors will also propagate to symbolic computations
because, for example, Equal relies on the significance arithmetic:

In[3]:=
Derivative[0, 0, 1, 0][Hypergeometric2F1][1, 1, 1, 1 - E] == 1/E

N::meprec: Internal precision limit $MaxExtraPrecision = 50.
  reached while evaluating -1/E +
  Derivative[0, 0, 1, 0][Hypergeometric2F1][1, 1, 1, 1 - E].

Out[3]=
False

In fact those two values are exactly equal; increasing
$MaxExtraPrecision to 100 gives the same incorrect answer even without
generating N::meprec message. An example involving MeijerG:

In[4]:=
N[MeijerG[{{0, 1}, {}}, {{-1/2, 1/2}, {}}, 1], 20]

<N::meprec warning>

Out[4]=
-9.8681345642064808273`5.348145572317817 + 0``4.353910509317925*I

The exact value is -Pi^2, which means that the last digit in the
output is incorrect; besides, setting $MaxExtraPrecision=1000 and
reevaluating this expression we obtain the same five significant
digits -- no increase in precision at all. Incidentally, here's the
result of FunctionExpand:

In[5]:=
MeijerG[{{0, 1}, {}}, {{-1/2, 1/2}, {}}, 1] // FunctionExpand

Infinity::indet: Indeterminate expression
  0*Pi^2*ComplexInfinity encountered.

Out[5]=
Indeterminate

FunctionExpand seems to give incorrect results for the hypergeometric
function and MeijerG rather often; here are a couple more examples:

In[6]:=
Hypergeometric2F1[1, 2, 1/2, -1/2] // FunctionExpand

Out[6]=
((6 + 3*I) - 2*Sqrt[3]*ArcCsch[Sqrt[2]])/9

(the correct value is (3-2*Sqrt[3]*ArcCsch[Sqrt[2]])/9).

In[7]:=
MeijerG[{{}, {1, 1}}, {{0, 0, 2}, {}}, 1] // FunctionExpand

Out[7]=
-1 + EulerGamma - (1 - E + E*EulerGamma - E*ExpIntegralEi[-1])/E

(here the sign is reversed).

Another example:

In[8]:=
N[MeijerG[{{-1, -1}, {}}, {{-3/2, -1/2}, {}}, 1], 20]

<N::meprec warning>

Out[8]=
3.4670435591329474158`5.506706408852548 + 0``4.966747110786506*I

This one is actually quite interesting, because it gives us an example
where MeijerG cannot be evaluated as the sum of residues; here's the
correct numerical value (the exact value is Pi^2/4):

In[9]:=
F[s_] = Gamma[2 - s]^2*Gamma[s - 3/2]*Gamma[s - 1/2];
1/(2*Pi*I)*NIntegrate[F[s], {s, -I*Infinity, 7/4, I*Infinity}]

Out[10]=
2.4674011002839356 + 0.*I

If we specify the contour of integration as above, the integration
will be performed along the vertical line Re[s]==7/4. And here's the
sum of the residues of F[s]:

In[11]:=
<<numericalmath`
Sum[NResidue[F[s], {s, k + 1/2}], {k, -100, 1}]

Out[12]=
3.4649411066986744 + 6.036664224051691*^-14*I

Just as well we could take the sum of the residues in the right
half-plane, then the result would be wrong by -1. The reason of course
is that it's simply mathematically incorrect, since as Abs[s] goes to
Infinity, F[s] does not tend to zero uniformly in Arg[s] and we cannot
neglect the integral along the big circle:

In[13]:=
1/(2*Pi*I)*NIntegrate[
  Evaluate[(F[s] /. s -> 7/4 + 1000*E^(I*t))*1000*I*E^(I*t)],
  {t, Pi/2, Pi - 10^-3, Pi + 10^-3, 3*Pi/2}]

Out[13]=
0.9997312041116979 - 2.9685177026992824*^-15*I

This is the missing unity; apparently, in this case Mathematica uses
the residue theorem in some form, hence the error.

Maxim Rytin
m.r at inbox.ru


  • Prev by Date: Re: what kind of a programming language is Mathematica?
  • Next by Date: Re: question about Integrate
  • Previous by thread: Re: customizing Integrate with Unprotect
  • Next by thread: neuroscience