Why such a big difference between results of Integrate and NIntegrate?

Studying an interesting article of Daniel Lichtblau, I consider a variation of an example from it, calculating an improper integral

Integrate[RealAbs[Sin[x - y]]^(-2/3), {x, 0, Pi}, {y, 0, Pi}]


-((12 Sqrt[\[Pi]] Gamma[-(1/3)] HypergeometricPFQ[{1/6, 1/2, 2/3}, {7/6, 7/6}, 1])/ Gamma[1/6])

This is not very useful analytic expression, so

N[%]


16.7126

Then I compare that result with the numeric one

NIntegrate[RealAbs[Sin[x - y]]^(-2/3), {x, 0, Pi}, {y, 0, Pi},
Exclusions -> {y == x}, AccuracyGoal -> 3, PrecisionGoal -> 3]


22.8915

The latest result is produced without any warning. How to explain so big difference between the numbers? Could this difference be decreased? Which result is more reliable?

• It should be fairly simple to reduce this to a 1D integral (x-y -> z). That ought to make it more obvious which is the correct answer. Apr 20, 2021 at 19:42
• @mikado: Thank you. However, this reduction is not so simple. Apr 20, 2021 at 19:45
• Note that you have exactly one period of the function over the integration interval. This makes the change of variable much simpler. In fact, you can just replace x-y with x Apr 20, 2021 at 19:47
• @mikado: Can you kindly ground your claim? How about the difference between the exact and numeric double integrals? Apr 20, 2021 at 19:50
• @mikado: Did you pay your attention that the set $x-y=0$ is smaller than the set $x-y=\pm \frac 1 2$ in rhe square $[0,2\pi]^2$? Apr 20, 2021 at 19:59

For a real function we can transform RealAbs[Sin[x - y]]^(-2/3) to 1/(1 - Cos[x - y]^2)^(1/3), then we have

Integrate[1/(1 - Cos[x - y]^2)^(1/3), {x, 0, Pi}, {y, 0, Pi}]

Out[]= (\[Pi]^(3/2) Gamma[1/6])/Gamma[2/3]


Numerical result is 22.88949310061915. It can be compare to

NIntegrate[1/(1 - Cos[x - y]^2)^(1/3), {x, 0, Pi}, {y, 0, Pi},
Exclusions -> {x == y}, WorkingPrecision -> 20]

Out[]= 22.889493100619489917

• MA11.1 gives your analytical result $\frac{\pi ^{3/2} \Gamma \left(\frac{1}{6}\right)}{\Gamma \left(\frac{2}{3}\right)}$ directly, without any reformulations. Apr 20, 2021 at 21:03
• @AlexTrounev: All that is not so simple. Let us consider the inner integral Integrate[1/(1 - Cos[x - y]^2)^(1/3), {y, 0, Pi}, Assumptions -> x >= 0 && x <= Pi] which results in ConditionalExpression[(6 Sqrt[\[Pi]] Gamma[7/6])/Gamma[2/3], 2 x < \[Pi]]. Pay your attention to 2 x < \[Pi]]. I am using 12.2 on Windows 10 Pro. Thank you anyway. Apr 20, 2021 at 21:10
• The command Integrate[1/(1 - Cos[x - y]^2)^(1/3), {y, 0, Pi}, Assumptions -> x >= Pi/2 && x <= Pi] results in (Sqrt[\[Pi]] Gamma[1/6])/Gamma[2/3], Combining these, we obtain the right (I hope) answer. Apr 20, 2021 at 21:22
• @AlexTrounev: Don't make me laudh. It's still unclear to me whether Integrate[1/(1 - Cos[x - y]^2)^(1/3), {y, 0, Pi}, Assumptions -> x >= Pi/2 && x <= Pi] formally applies the Newton-Leibniz formula. not paying attention to the singularity at y=x or the command does consider it as an improper integral. NIntegrate also may switch to the symbolic integration. BTW, Integrate[1/(1 - Cos[x - y]^2)^(1/3), y, Assumptions -> x >= Pi/2 && x <= Pi results in -((3 Sqrt[Cos[x - y]^2] Hypergeometric2F1[1/6, 1/2, 7/6, Sin[x - y]^2] Tan[x - y])/(Sin[ x - y]^2)^(1/3)) . Regard. Apr 20, 2021 at 21:37
• @AlexTrounev: The question is still open. Do you understand me? Don't hesitate to ask for further explanation in need. Regard again. Apr 20, 2021 at 21:40

Here is my answer to my question. First, it is enough to consider the integration over the triangle with vertices at $$(0,0),\,(0,\pi),\,(\pi,\pi)$$. Second, we make the change of the variables $$s=y-x,\,t=x+y$$, using the idea of @mikado. The inverse map is $$x =\frac {t-s} 2,\, y=\frac {t+s} 2$$ and its Jacobian determinant equals $$-\frac 1 2$$. Now the set of the integration is the triangle with its vertices at $$(0,0),\,(0,2\pi),\,(\pi,\pi)$$ and the double integral under consideration becomes $$\int\limits_0^\pi \int\limits_s^{2\pi-s} \frac 1 {\sin(s)^{2/3}}\left|- \frac 1 2 \right| \,dtds= \int\limits_0^\pi \frac {\pi-s} {\sin(s)^{2/3}}\,ds.$$

Mathematica calculates the latter integral by

Integrate[(Pi - s)/Sin[s]^(2/3), {s, 0, Pi}]


(3 \[Pi]^(3/2) Gamma[7/6])/Gamma[2/3]

N[(3 \[Pi]^(3/2) Gamma[7/6])/Gamma[2/3]]


11.4447

This coinsides with the result by @AlexTrounev. Let us consider how Mathematica 12.2 does the job. Its result for the antiderivative

Integrate[(Pi - s)/Sin[s]^(2/3), s, Assumptions -> s > 0 && s <= Pi]


1/2 Sin[s]^( 1/3) (-6 s Cos[s] Hypergeometric2F1[2/3, 1, 7/6, Sin[s]^2] + (Sqrt[\[Pi]] Gamma[1/3] HypergeometricPFQ[{2/3, 2/3, 1}, {7/6, 5/3}, Sin[s]^2] Sin[s])/(2^(1/3) Gamma[7/6] Gamma[5/3]) - (2 \[Pi] Cos[s] Hypergeometric2F1[1/2, 5/6, 3/2, Cos[s]^2])/(Sin[s]^2)^(1/6))

is not correct in view of

Plot[%,{s,0,Pi}]


which shows a discontinuous function.

I think Mathematica knows the integral Integrate[(Pi - s)/Sin[s]^(2/3), {s, 0, Pi}] as a table value, though I don't find it in a handbook of Gradshtein&Ruezhik. Numeric calculations confirm it by

Table[NIntegrate[(Pi - s)/Sin[s]^(2/3), {s, 10^(-k), Pi}], {k, 1, 4}]


{7.10428, 9.41585, 10.5023, 11.0073}

• Do you really think that Mathematica knows? There are rules for interpretation every combination of parameters. The question is why these interpretations are differ for Integrate and NIntegrate in versions 12.0-12.2 while in v.11.1 they are the same? My opinion that it is connected with interpretation of RealAbs. It is why I remove this function in my answer. The second problem is interpretation of power 2/3. Apr 21, 2021 at 11:17
• @AlexTrounev: Your guess that RealAbs causes an incorrect answer is built on the sand. The command Integrate[RealAbs[Sin[x - y]]^(-1/2), {x, 0, Pi}, {y, 0, Pi}] produces the correct answer. Apr 21, 2021 at 11:33
• @AlexTrounev: I'd like to add that Integrate[(Pi - s)/RealAbs[Sin[s]]^(2/3), {s, 0, Pi}] also produces (3 \[Pi]^(3/2) Gamma[7/6])/Gamma[2/3]. Apr 21, 2021 at 11:43
• Even your analysis is good (+1), the main problem is not explained. Why Integrate and NIntegrate give different answer? Apr 21, 2021 at 12:59
• @AlexTrounev: After the change of variables or/and your (and MichaelE2's) representation of the integrand the results are concordant. I don't know " How to explain so big difference between the numbers?" so I asked it. Apr 21, 2021 at 15:56

I'm left unsure what is the correct answer from the above analyses. Personally, I believe the integral should be evaluated as an improper or principal-valued integral. In this case, the singularity is the line y=x. Below I create a principal-valued integral function f1[x], then integrate this function over the requisite domain:

f1[x_?NumericQ] :=
NIntegrate[RealAbs[Sin[x - y]]^(-2/3), {y, 0, x,
Pi},Method -> "PrincipalValue"]
NIntegrate[f1[x], {x, 0, Pi}]

Out[8]= 22.8895


It's just more simple-understanding this way and one more intuitively reasonable I think. Also I have doubt change of variable is applicable in a domain with singularities.

• Why Integrate and NIntegrate compute different output? Apr 21, 2021 at 13:07
• I reviewed briefly Daniel's article. There can be problems with Integrate handling singularities. I suspect it is not handling them correctly in this case. There is also the problem of interpreting associated multi-valued functions involved. That jump in the plot above is a clear example of a branch-discontinuity which would lead to an incorrect answer if an analytically-continuous result is desired. Apr 21, 2021 at 13:20
• It is very good. But why in your answer there is no example with Integrate? Apr 21, 2021 at 13:24
• Good point Alex. However, Integrate returns an analytic expression containing the underlying HypergeometricPFQ function which is multivalued. In order to use this function correctly, an analytically-continuous route over the branching would be needed similar to evaluating $$\oint 1/z dz$$ in terms of the logarithmic antiderivative. Not sure how to do this for Hypergeometric and also I think Daniel and many others here more qualified to analyze this than me. Apr 21, 2021 at 13:32
• The double improper integral under consideration converges, so Method->"PricipalValue" does nothing. The same result is produced by f1[x_?NumericQ] := NIntegrate[RealAbs[Sin[x - y]]^(-2/3), {y, 0, Pi}] NIntegrate[f1[x], {x, 0, Pi}]. You simply numerically calculate the improper integral in a complicated way, no more and no less. This is not it. Apr 21, 2021 at 13:39