2
$\begingroup$

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?

$\endgroup$

1 Answer 1

4
$\begingroup$

Your first integral might be solved in two steps!

As far as I know, Finite Element Meshing is only provided up to dimension 3 in Mathematica. The error message "...method FiniteElement " led me to a stepwise approach, which works without error message:

int[p_?NumericQ, q_?NumericQ, r_?NumericQ] := 
NIntegrate[1/Sqrt[(x - p)^2 + (y - q)^2 + (z - r)^2], {x, y, z} \[Element]Ball[{0, 0, 0}, 1] ]    

NIntegrate[int[p , q , r ], {p, q, r} \[Element] Ball[{0, 0, 0}, 1]]
(*21.0552*)
$\endgroup$
2
  • $\begingroup$ My first feedback: your code is running without any result on my comp during more than15 minutes. $\endgroup$
    – user64494
    Apr 27, 2021 at 11:26
  • $\begingroup$ +1. Thank you: it works. Few remarks.(i) The execution takes approximately 30 minutes on my comp. (ii) It does not work for the four-dimensional unit ball. (iii) Don't understand too much difference between 21.0552 produced by your code and 20.32729987338035891791629 produced by AdaptiveMonteCarl with AccuracyGoal -> 2, PrecisionGoal -> 2, WorkingPrecision -> 25. If I correctly understand it, this option guarrantees the relative error 0.01. $\endgroup$
    – user64494
    Apr 27, 2021 at 11:58

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.