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},
Method -> "AdaptiveMonteCarlo"]
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?
