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.12327556039856680489716`25. and 0.2251832352724191939747851`

25. 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`

witout it.

The questions arise: what are other methods to calculate such integrals?

is a three-digit result possible with MonteCarlo methods?