# Problems with improper integrals in higher dimensions

Calculating numerically certain energy integrals in three and four dimensions related to the Riesz potential and the capacity, I try

NIntegrate[ 1/Sqrt[(x - p)^2 + (y - q)^2 + (z - r)^2], {x, y, z} \[Element]
Ball[{0, 0, 0}, 1], {p, q, r} \[Element] Ball[{0, 0, 0}, 1],
AccuracyGoal -> 3, PrecisionGoal -> 3, WorkingPrecision -> 25,
Exclusions -> {(x - p)^2 + (y - q)^2 + (z - r)^2 == 0}]

0

and a warning message "NIntegrate::moptxn: The option SymbolicProcessing of the method FiniteElement is not one of {Method,MeshOptions}." As far as I understand it, this means that the NIntegrate command does not accept the set of the integration in the form {x, y, z} \[Element] Ball[{0, 0, 0}, 1], {p, q, r} \[Element] Ball[{0, 0, 0}, 1]. However, if it is so, the input, not 0, should be returned.

My next try is

NIntegrate[ 1/((x - p)^2 + (y - q)^2 + (z - r)^2)^(1/2), {x, -1, 1}, {y, -Sqrt[1 - x^2], Sqrt[1 - x^2]},
{z, -Sqrt[1 - x^2 - y^2], Sqrt[1 - x^2 - y^2]}, {p, -1, 1}, {q, -Sqrt[1 - p^2], Sqrt[1 - p^2]},
{r, -Sqrt[1 - p^2 - q^2], Sqrt[1 - p^2 - q^2]}, AccuracyGoal -> 3, PrecisionGoal -> 3, WorkingPrecision -> 25,
Exclusions -> {(x - p)^2 + (y - q)^2 + (z - r)^2 == 0}]

-1244.482640558337558417913

and a warning "NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 54 recursive bisections in z near {x,y,z,p,q,r} = <<1>>. NIntegrate obtained -1244.482640558337558417913 and 3534.518334552443660338233`25. for the integral and error estimates." A similar issue with

NIntegrate[1/((x - p)^2 + (y - q)^2 + (z - r)^2)^(1/2)*
Boole[x^2 + y^2 + z^2 <= 1 && p^2 + q^2 + r^2 <= 1], {x, -1, 1}, {y, -1, 1}, {z, -1, 1}, {p, -1, 1},
{q, -1, 1}, {r, -1, 1}, AccuracyGoal -> 3, PrecisionGoal -> 3, WorkingPrecision -> 25,
Exclusions -> {(x - p)^2 + (y - q)^2 + (z - r)^2 == 0}]

-1203.034524853306755966531

I wonder negative numbers since the integrand is positive. Then i switch to MonteCarlo methods.

NIntegrate[ 1/((x - p)^2 + (y - q)^2 + (z - r)^2)^(1/2)*
Boole[x^2 + y^2 + z^2 <= 1 && p^2 + q^2 + r^2 <= 1], {x, -1, 1}, {y, -1, 1}, {z, -1, 1}, {p, -1, 1},
{q, -1, 1}, {r, -1, 1}, AccuracyGoal -> 2, PrecisionGoal -> 2, WorkingPrecision -> 25,
Exclusions -> {(x - p)^2 + (y - q)^2 + (z - r)^2 == 0},Method -> "QuasiMonteCarlo"]

21.12327556039856680489716

and a warning "NIntegrate::maxp: The integral failed to converge after 50000 integrand evaluations. NIntegrate obtained 21.1232755603985668048971625. and 0.225183235272419193974785125. for the integral and error estimates.". A more or less reliable result is obtained by

NIntegrate[ 1/((x - p)^2 + (y - q)^2 + (z - r)^2)^(1/2)*
Boole[x^2 + y^2 + z^2 <= 1 && p^2 + q^2 + r^2 <= 1], {x, -1, 1}, {y, -1, 1}, {z, -1, 1}, {p, -1, 1},
{q, -1, 1}, {r, -1, 1},  AccuracyGoal -> 2, PrecisionGoal -> 2, WorkingPrecision -> 25,
Exclusions -> {(x - p)^2 + (y - q)^2 + (z - r)^2 == 0}, Method -> "AdaptiveMonteCarlo"]

20.32729987338035891791629

without any warning. Unfortunately, AccuracyGoal -> 3, PrecisionGoal -> 3 is not achieved. Also this works in eight dimensions:

NIntegrate[1/((x - p)^2 + (y - q)^2 + (z - r)^2 + (w - s)^2)^((4 - 2)/2)*
Boole[x^2 + y^2 + z^2 + w^2 <= 1 &&  p^2 + q^2 + r^2 + s^2 <= 1], {x, -1, 1}, {y, -1, 1}, {z,-1, 1},
{w, -1, 1}, {p, -1, 1}, {q, -1, 1}, {r, -1, 1}, {s, -1, 1}, AccuracyGoal -> 2, PrecisionGoal -> 2,
WorkingPrecision -> 25,Exclusions -> {(x - p)^2 + (y - q)^2 + (z - r)^2 + (w - s)^2 == 0},

25.78510573458365881830859

BTW, the Exclusions option works in the above: compare with 24.68762920929857902438239 without it.

The questions arise: what are other methods to calculate such integrals? is a three-digit result possible with MonteCarlo methods?