The question is about Mathematica's result for the integral
$$ \int_0^\infty \frac{16}{\pi^2} \left( \frac{\sqrt{\pi}}{4} r^2 % G_{1,3}^{2,1}\left(\frac{r^2}{4}\Bigg| \begin{array}{l} 0 \\ 0,0,-\frac{1}{2} \\ \end{array} \right) % \right)^2 dr \,,$$
whereby the integrand will be denoted by $f(r)$ in the following. $G$ denotes the Meijer G-function (name of the corresponding Mathematica function: MeijerG). According to Mathematica's Integrate function this integral is $1$. However, with NIntegrate one obtains the following results:
$$ \int_0^5 f(r) d r > 7 \,, $$ $$\int_0^{10} f(r) d r > 16 \,. $$ Furthermore Mathematica confirms that $$ f(r) \geq 0 \,,\ \ \forall r \geq 0 $$ holds. Using this together with the results from the numeric integrals one obtains $$ \int_0^\infty f(r) d r = \int_0^{10} f(r) d r + \int_{10}^\infty f(r) d r > 16 + \int_{10}^\infty f(r) d r > 16 \,,$$ since the second integral has to be $\geq 0$ due to the non-negativity of the integrand. This contradicts the result from Mathematica's Integrate, which says that the integral is $1$.

The question is: Where is the fault leading to the contradiction? Which statement is true?
(I would naively assume that Integrate is wrong in this case.)

Note: In the case one doubts the results for the numeric integrals from NIntegrate, one can plot the integrand. These plots make the results seem very plausible. A screenshot of the plots and of the used Mathematica commands can be found below. Additionally, the used code is listed below for the case that some wants to test it.

Screenshot of the Mathematica notebook where the integrals are evaluated

Used Mathematica code

intgd[r_] := 16/\[Pi]^2*(1/4 Sqrt[\[Pi]] r^2 MeijerG[{{0}, {}}, {{0, 0}, {-(1/2)}}, r^2/4])^2
Integrate[intgd[r], {r, 0, Infinity}]
NIntegrate[intgd[r], {r, 0, 5}]
FullSimplify[intgd[r] >= 0, Assumptions -> {r >= 0}]
GraphicsRow[{Plot[intgd[r], {r, 0, 20}], Plot[intgd[r], {r, 0, 50}]}]

The used Mathematica version is 12.2.


1 Answer 1



(* "12.3.1 for Mac OS X x86 (64-bit) (June 19, 2021)" *)

intgd[r_] := 
 16/π^2*(1/4 Sqrt[π] r^2 MeijerG[{{0}, {}}, {{0, 0}, {-(1/2)}}, 

lim = Limit[intgd[r], r -> Infinity]

(* 16/π^2 *)

For larger values of r, arbitrary-precision rather than machine precision is required for the Plot

Plot[intgd[r], {r, 0, 50},
 WorkingPrecision -> 30,
 PlotRange -> All,
 GridLines -> {None, {lim}},
 GridLinesStyle -> Directive[Red, AbsoluteThickness[1]]]

enter image description here

An alternate representation of the integrand for r > 0 is

intgd2[r_] = intgd[r] // FunctionExpand // Simplify[#, r > 0] &

(* (16 r^2 (CoshIntegral[r] Sinh[r] - Cosh[r] SinhIntegral[r])^2)/π^2 *)

intgd[r] == intgd2[r] // FullSimplify[#, r > 0] &

(* True *)

For this alternate representation, even greater precision is required in the Plot

Plot[intgd2[r], {r, 0, 50},
 WorkingPrecision -> 50,
 PlotRange -> All,
 GridLines -> {None, {lim}},
 GridLinesStyle -> Directive[Red, AbsoluteThickness[1]]]

enter image description here

As expected, the integral on {0, Infinity} does not converge

Integrate[intgd2[r], {r, 0, Infinity}]

enter image description here

  • $\begingroup$ So, this is a nice way to show that Integrate[intgd[r], {r, 0, Infinity}] produces an incorrect result. This implies that Mathematica's Integrate has a bug here, right? (And thanks for showing how to plot the integrand adequately!) $\endgroup$ Commented Nov 5, 2021 at 20:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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