4
$\begingroup$

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:

enter image description here

How can I get Mathematica to solve this problem numerically?

$\endgroup$
  • 1
    $\begingroup$ 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? $\endgroup$ – John Doty Mar 25 at 12:20
  • $\begingroup$ 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. $\endgroup$ – Javier Alaminos Mar 25 at 12:31
  • 2
    $\begingroup$ Use option Exclusions -> {-1, 1, y + x == 0}] $\endgroup$ – user18792 Mar 25 at 12:43
  • 1
    $\begingroup$ All that is built on the sand because the PrincipalValue option for multivariate integrals is undocumented. $\endgroup$ – user64494 Mar 25 at 20:56
2
$\begingroup$

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

enter image description here

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

Have fun!

$\endgroup$
  • 3
    $\begingroup$ I don't understand what you're plotting. The value of this integral is 0 and your result is around 12. $\endgroup$ – Javier Alaminos Mar 25 at 14:10
  • $\begingroup$ You are right, it is not the same integral, since I took Sqrt[(x-y)^2+eps^2] instead of x-y. $\endgroup$ – Alexei Boulbitch Mar 25 at 15:52
  • $\begingroup$ So, to solve the original integral what do I have to do? $\endgroup$ – Javier Alaminos Mar 25 at 18:28
  • $\begingroup$ 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. $\endgroup$ – Michael Seifert Mar 25 at 18:35
2
$\begingroup$

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.}, 
  Method -> "LocalAdaptive"]

=>

7.73766*10^-17

For now, we can see that the integral value is close to zero.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Javier Alaminos Mar 25 at 18:05
  • $\begingroup$ @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. $\endgroup$ – Xminer Mar 25 at 20:02

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.