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