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