3
$\begingroup$

I was asked to look at a complicated definite integral that could be integrated analytically, but numerically did not behave as expected. After a lot of simplifications I arrived at the following:

In[4]:= int= Integrate[x BesselJ[0, 1/x], {x, 0, 1/2000}]

Out[4]= -(1999999/8000000)+EulerGamma/4-
        31250 HypergeometricPFQ[{1,1},{2,3,3},-1000000]+(3 Log[10])/4

When we inspect the plot of the integrand, we see that the absolute value of the integral must be considerably less than $2\times 10^{-5}$. We find:

In[5]:= N[int]

Out[5]= -2.04448*10^-12

In[6]:= N[int, 5]

Out[6]= 1.6212

The second result seems to be wrong. The culprit is the hypergeometric function:

In[7]:= N[ HypergeometricPFQ[{1,1},{2,3,3},-1000000]]

Out[7]= 0.0000518798

In[8]:= N[ HypergeometricPFQ[{1,1},{2,3,3},-1000000], 5]

Out[8]= 6.5424*10^-17

I also observed the very frustrating behaviour that I cannot always reproduce the incorrect output shown above. At the moment, Mathematica behaves stable: when I evaluate the above cells once again, I obtain the same output. But not rarely I observed that the results of N with a second argument suddenly were in agreement with the output without a second argument. To see the erroneous behavior again, I had to change the upper limit in the integral a little bit. I am running version 11.3 under Windows 10.

Am I right that the above arbitrary precision result for HypergeometricPFQ is incorrect? Can somebody reproduce the incorrect result?

$\endgroup$
  • $\begingroup$ NIntegrate[x BesselJ[0, 1/x], {x, 0, 1/2000}, WorkingPrecision -> 32] suggests the arb. prec. result is wrong. and the machine precision result is accurate. (V11.3.0) $\endgroup$ – Michael E2 Jul 2 '18 at 18:00
  • $\begingroup$ Also suggested by N[HypergeometricPFQ[{1, 1}, {2, 3, 3}, -1000000], 100] $\endgroup$ – ilian Jul 2 '18 at 19:23
  • $\begingroup$ @Ilian. Your result shows that you arrived at the correct result. My frustration is that sometimes, I do not know when, I get this correct result as well. Did you see the incorrect result? $\endgroup$ – Fred Simons Jul 3 '18 at 7:06
  • $\begingroup$ What version are you using? In Version 11.3 I get consistent and correct answers for all these cases. $\endgroup$ – TheDoctor Jul 8 '18 at 2:31
  • $\begingroup$ @TheDoctor, Thanks. It looks like I am the only one getting this incorrect result (V11.3, windows). I reported it to WRI. $\endgroup$ – Fred Simons Jul 8 '18 at 8:09
3
$\begingroup$

Here is a workaround. The integral under consideration, after making a substitution, can be expressed as

$$\int_{2000}^\infty\frac{J_0(u)}{u^3}\mathrm du$$

Through Mellin transform techniques, this can be converted to a Meijer $G$ expression:

$$\frac{1}{8}\operatorname{G}_{1,3}^{2,0}\left(\frac{2000}{2},\frac{1}{2}\middle |{{1}\atop{-1,0,-1}}\right)$$

or in Mathematica syntax, MeijerG[{{}, {1}}, {{-1, 0}, {-1}}, 2000/2, 1/2]/8. Using FunctionExpand[] on it confirms the closed form given in the OP.

This seems to work with arbitrary precision, unlike the hypergeometric expression:

N[MeijerG[{{}, {1}}, {{-1, 0}, {-1}}, 2000/2, 1/2]/8, 5]
   -2.0445*10^-12

Alternatively, a slight change in the Mellin route yields the equivalent Meijer $G$ expression

$$\frac{1}{8}\operatorname{G}_{3,1}^{0,2}\left(\frac{2}{2000},\frac{1}{2}\middle |{{1,2,2}\atop{0}}\right)$$

(in Mathematica form, MeijerG[{{1, 2}, {2}}, {{}, {0}}, 2/2000, 1/2]/8) which is related to the first one via the "flip identity". This, too, seems to also work for arbitrary precision evaluations.

$\endgroup$
  • $\begingroup$ Just back from a short holiday, I found your answer. That is a wonderful piece of mathematics! I forwarded your answer to the person who asked me to look at the integral. $\endgroup$ – Fred Simons Oct 1 '18 at 16:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.