3
$\begingroup$

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?

Improve this question
$\endgroup$
12
  • $\begingroup$ 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. $\endgroup$
    – mikado
    Apr 20, 2021 at 19:42
  • $\begingroup$ @mikado: Thank you. However, this reduction is not so simple. $\endgroup$
    – user64494
    Apr 20, 2021 at 19:45
  • $\begingroup$ 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 $\endgroup$
    – mikado
    Apr 20, 2021 at 19:47
  • $\begingroup$ @mikado: Can you kindly ground your claim? How about the difference between the exact and numeric double integrals? $\endgroup$
    – user64494
    Apr 20, 2021 at 19:50
  • $\begingroup$ @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$? $\endgroup$
    – user64494
    Apr 20, 2021 at 19:59

3 Answers 3

Reset to default
3
$\begingroup$

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
Improve this answer
$\endgroup$
8
  • 1
    $\begingroup$ 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. $\endgroup$
    – yarchik
    Apr 20, 2021 at 21:03
  • $\begingroup$ @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. $\endgroup$
    – user64494
    Apr 20, 2021 at 21:10
  • $\begingroup$ 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. $\endgroup$
    – user64494
    Apr 20, 2021 at 21:22
  • $\begingroup$ @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. $\endgroup$
    – user64494
    Apr 20, 2021 at 21:37
  • $\begingroup$ @AlexTrounev: The question is still open. Do you understand me? Don't hesitate to ask for further explanation in need. Regard again. $\endgroup$
    – user64494
    Apr 20, 2021 at 21:40
1
$\begingroup$

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}]

enter image description here 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}

Improve this answer
$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$
    – Alex Trounev
    Apr 21, 2021 at 11:17
  • $\begingroup$ @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. $\endgroup$
    – user64494
    Apr 21, 2021 at 11:33
  • $\begingroup$ @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]. $\endgroup$
    – user64494
    Apr 21, 2021 at 11:43
  • $\begingroup$ Even your analysis is good (+1), the main problem is not explained. Why Integrate and NIntegrate give different answer? $\endgroup$
    – Alex Trounev
    Apr 21, 2021 at 12:59
  • $\begingroup$ @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. $\endgroup$
    – user64494
    Apr 21, 2021 at 15:56
1
$\begingroup$

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.

Improve this answer
$\endgroup$
14
  • $\begingroup$ Why Integrate and NIntegrate compute different output? $\endgroup$
    – Alex Trounev
    Apr 21, 2021 at 13:07
  • 1
    $\begingroup$ 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. $\endgroup$
    – Dominic
    Apr 21, 2021 at 13:20
  • $\begingroup$ It is very good. But why in your answer there is no example with Integrate? $\endgroup$
    – Alex Trounev
    Apr 21, 2021 at 13:24
  • $\begingroup$ 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. $\endgroup$
    – Dominic
    Apr 21, 2021 at 13:32
  • $\begingroup$ 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. $\endgroup$
    – user64494
    Apr 21, 2021 at 13:39

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.