# NIntegrate: How can I solve this integral numerically? NIntegrate fails while Integrate works

I know the exact solution of the principal value of this integral is equal to zero:

$$\int_{-1}^{1}\int_{-1}^{1}\frac{x^2}{\sqrt{1-x^2}}\frac{\sqrt{1-y^2}}{y-x}dydx=0$$

doing:

Integrate[x^2/Sqrt[1 - x^2] Sqrt[1 - y^2]/(y - x), {x, -1, 1}, {y, -1, 1},
PrincipalValue -> True]


but I want to do it numerically and it doesn't work:

NIntegrate[x^2/Sqrt[1 - x^2] Sqrt[1 - y^2]/(y - x), {x, -1, 1}, {y, -1, 1}]


This is the error message returned:

How can I get Mathematica to solve this problem numerically?

• The issue is that the integrand approaches infinity as x->±1, x->y, and y->x. That kind of behavior is toxic to numerical methods: you need to reason out a way to deal with it, not merely probe it numerically. PrincipalValue -> True gives you access to automated reasoning in this case, and you've solved the problem that way. Do you have a different problem you're trying to solve? – John Doty Mar 25 at 12:20
• Yes, I'm trying to solve a similar integral, when x^2 is multplied by exp^(-i*b*(x + y)). So, firstly I wanted to try to solve this known integral. – Javier Alaminos Mar 25 at 12:31
• Use option Exclusions -> {-1, 1, y + x == 0}] – user18792 Mar 25 at 12:43
• All that is built on the sand because the PrincipalValue option for multivariate integrals is undocumented. – user64494 Mar 25 at 20:56

The main problem is the point x=y. In principle, it seems that there the integral is singular. If you agree to get the principal value of it, you may exclude this point by a regularization as follows

NIntegrate[
x^2/Sqrt[1 - x^2] Sqrt[1 - y^2]/Sqrt[(y - x)^2 + i^-2], {x, -1 + 1/i,
1 - 1/i}, {y, -1 + 1/i, 1 - 1/i}, Method -> "AdaptiveMonteCarlo"]


where i is a large number. Then you may increase i and check the convergence of the integral:

    lst = Table[{1/i,
NIntegrate[
x^2/Sqrt[1 - x^2] Sqrt[1 - y^2]/
Sqrt[(y - x)^2 + i^-2], {x, -1 + 1/i, 1 - 1/i}, {y, -1 + 1/i,
1 - 1/i}, Method -> "AdaptiveMonteCarlo"]}, {i, 1000, 100000,
1000}] // N;

ListLogPlot[lst /. {x_, y_} -> {1/x, y}, Frame -> True,
FrameLabel -> {Style["Number i", 16], Style["Integral", 16]}]


yielding this

One can further a few other methods which may eventually enable a more accurate estimate of the integral.

Have fun!

• I don't understand what you're plotting. The value of this integral is 0 and your result is around 12. – Javier Alaminos Mar 25 at 14:10
• You are right, it is not the same integral, since I took Sqrt[(x-y)^2+eps^2] instead of x-y. – Alexei Boulbitch Mar 25 at 15:52
• So, to solve the original integral what do I have to do? – Javier Alaminos Mar 25 at 18:28
• Note that $\sqrt{(y-x)^2 + \epsilon^2}$ approaches $|y - x|$ as $\epsilon \to 0$, not the $y - x$ that's in the original integrand. – Michael Seifert Mar 25 at 18:35

As the others say,simply integrate by avoiding singular points?

Fixed.

# Try other integral.

target = Compile[{{x, _Real}, {y, _Real}},
x/\[Sqrt](1 - x) \[Sqrt](1 - y^2)/(y - x)];
integrand[x_, y_] := If[Or[(1 - x) == 0, y == x], 0, target[x, y]];
Quiet@NIntegrate[integrand[x, y], {x, -1., 1.}, {y, -1., 1.}]


=>

-4.06259


Integration by manual

.

Plus @@ Flatten@
Table[integrand[x, y]*0.001*0.001, {x, -1., 1., 0.001}, {y, -1., 1.,
0.001}]


=>

-3.99866


Integration by NIntegrate

N@Integrate[
x/Sqrt[1 - x] Sqrt[1 - y^2]/(y - x), {x, -1, 1}, {y, -1, 1},
PrincipalValue -> True]


=>

-4.14669


# the question's integral.

target = Compile[{{x, _Real}, {y, _Real}},
x^2/\[Sqrt](1 - x^2) \[Sqrt](1 - y^2)/(y - x)];
integrand[x_, y_] := If[Or[(1 - x^2) == 0, y == x], 0, target[x, y]];
Quiet@NIntegrate[integrand[x, y], {x, -1., 1.}, {y, -1., 1.}]


=>

-0.4542


By manual.

Plus @@ Flatten@
Table[integrand[x, y]*0.1*0.1, {x, -1., 1., 0.1}, {y, -1., 1., 0.1}]


=>

-8.88178*10^-16


By other method.

Quiet@NIntegrate[integrand[x, y], {x, -1., 1.}, {y, -1., 1.},

7.73766*10^-17

• But I can't avoid the singularity $x==y$, because for example if I have $x$ instead of $x^2$ the result of integral is $-\pi^2/2$, and with your code the result is always 0. – Javier Alaminos Mar 25 at 18:05
• @JavierAlaminos　as your point,my code was always 0. it's mainly because my code returns Nothing when the condition met,I think. so just modified. – Xminer Mar 25 at 20:02