# Numerical integration of multidimensional integral

I need to compute quite a few multidimensional integrals with Bessel functions. Here is one example

NIntegrate[
Subscript[r, 1] Subscript[r, 2] Subscript[r, 3]
BesselK[0, Abs[Subscript[r, 1]]]^2 BesselK[0,
Abs[Subscript[r, 1] -
E^(I Subscript[\[Phi], 2]) Subscript[r, 2]]] BesselK[0,
Abs[E^(I Subscript[\[Phi], 2]) Subscript[r, 2] -
E^(I Subscript[\[Phi], 3]) Subscript[r,
3]]]^3 , {Subscript[\[Phi], 2], 0, 2 \[Pi]}, {Subscript[\[Phi],
3], 0, 2 \[Pi]}, {Subscript[r, 1], 0, \[Infinity]}, {Subscript[r,
2], 0, \[Infinity]}, {Subscript[r, 3], 0, \[Infinity]}]


So the accuracy of numerical integration is pretty bad with several of the methods I have tried. For example the standard GlobalAdaptive estimates an error about 15% and increasing MaxErrorIncreases from 2000 up to 10000 does not seem to help. I would be satisfied to get a correct answer within 1% accuracy and since I have to do many integrations of this type possibly in reasonable time.

The particular integral above in symbolic terms is (up to constants) $$\iiint d^2x_1 d^2x_2d^2x_3\,\, K_0(|x_1|)^2K_0(|x_1-x_2|)K_0(|x_2-x_3|)^3$$ The two-dimensional integrals all run over all 2d space with standard polar measure $$d^2x=r\,dr\,d\phi$$ and due to rotational symmetry integration over one of the angles can be eliminated. $$K_0(r)$$ is the modified Bessel function BesselK[0,r].

The integrand is everywhere positive so the integral is non-zero. Bessel function $$K_0(r)$$ has a logarithmic singularity at $$r=0$$ which is integrable in 1d and easily integrable in 2d since $$d^2x K_0(r) \simeq dr \,r\,K_0(r)$$. I suspect that it is this singularity that leads to very imprecise numerical integration, but I'm not sure what I can do about it.

## Edit 1

I'm starting a bounty because the accuracy of the method suggested by Ulrich Neumann does not seem to be sufficient and even well under control.

realabs = Sqrt[#^2] &
NIntegrate[r1 r2 r3 BesselK[0, realabs[r1]]^2 BesselK[0,Abs[r1 - E^(I \[Phi]2) r2]] BesselK[0,Abs[E^(I \[Phi]2) r2 - E^(I \[Phi]3) r3]]^3, {\[Phi]2, 0, 2 \[Pi]}, {\[Phi]3, 0, 2 \[Pi]}, {r1, 0, \[Infinity]}, {r2,0, \[Infinity]},{r3, 0, \[Infinity]}
, WorkingPrecision -> 15 , AccuracyGoal -> 5,Method -> {"AdaptiveMonteCarlo","SymbolicProcessing" -> 0}] // AbsoluteTiming
(*{60.5697, 7.37415114501064}*)


On different runs the result fluctuates within about 10%. What is worse is that the exact answer computed by the method suggested by Roman is something like

NIntegrate[(2 \[Pi])^2 r1 r2 r3 BesselK[0, r1]^2 BesselK[0,
r2] BesselK[0, r3]^3, {r1, 0, \[Infinity]}, {r2,
0, \[Infinity]}, {r3, 0, \[Infinity]}]
(* 11.5667 *)


which means that multidimensional numerical integration gives almost 50% error.

## Edit 2

I've tested suggestion by Alex Trounev on the following integral

g = Compile[{{r2, _Real}, {r3, _Real}, {r4, _Real}, {\[Phi]3, _Real}, \
{\[Phi]4, _Real}},
2 \[Pi] Abs[r2]^4 BesselK[0, r2] BesselK[0,
Abs[E^(I \[Phi]3) r3]]^2 BesselK[0,
Abs[r2 - E^(I \[Phi]4) r4]]^2 BesselK[0,
Abs[E^(I \[Phi]3) r3 - E^(I \[Phi]4) r4]] r2 r3 r4,
Parallelization -> True];
zero = $MachineEpsilon; L = Infinity; NIntegrate[ g[r2, r3, r4, \[Phi]3, \[Phi]4], {\[Phi]3, 0, 2 \[Pi]}, {\[Phi]4, 0, 2 \[Pi]}, {r2, zero, Infinity}, {r3, zero, Infinity}, {r4, zero, Infinity}, Method -> "QuasiMonteCarlo"]  Although using compiled function definetely speeds up the process the result does not seem to be robust. I do not know the exact value of this integral but the Monte-Carlo methods seem to fluctuate a lot. QuasiMonteCarlo gives about 10% error estimate. Is there a way to nail the value down more reliably? ## Edit 3 So the above integral seems to be captured well by the MonteCarlo methods. To verify this we could run NIntegrate[ g[r2, r3, r4, \[Phi]3, \[Phi]4], {\[Phi]3, 0, 2 \[Pi]}, {\[Phi]4, 0, 2 \[Pi]}, {r2, zero, Infinity}, {r3, zero, Infinity}, {r4, zero, Infinity}, Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 10000}]  which returns 84.14 with 0.35 error estimate, within $$1\%$$ as requested. ## 3 Answers This integral can be done analytically by substitution. Set $$x_1=r_1 e^{i\phi_1}$$, $$x_2=x_1+r_2 e^{i\phi_2}$$, and $$x_3=x_2+r_3 e^{i\phi_3}$$. With these definitions, $$|x_1|=r_1$$, $$|x_1-x_2|=r_2$$, and $$|x_2-x_3|=r_3$$: the integral separates into three, Integrate[(2π r1) (2π r2) (2π r3) r1^2 BesselK[0, r1]^2 BesselK[0, r2] BesselK[0, r3]^3, {r1, 0, ∞}, {r2, 0, ∞}, {r3, 0, ∞}] (* 4/3 π^(7/2) MeijerG[{{0, 0}, {1/2}}, {{0, 0, 0}, {}}, 4] *) % // N (* 48.4506 *)  Of course it's a stretch to call this result "analytic" when it contains a Meijer G function. However, it is now reduced to a standardized function, for which there are accurate and proven algorithms and there is no need for numerically optimized integration. Further, with FunctionExpand we can find the definition of this Meijer G function in terms of logarithms and polylogarithms: %% // FunctionExpand // FullSimplify (* 4/27 π^3 (Sqrt[3] π Log[4] - PolyGamma[1, 1/3] + PolyGamma[1, 2/3] + 6 I Sqrt[3] (PolyLog[2, 1/4 - (I Sqrt[3])/4] - PolyLog[2, 1/4 + (I Sqrt[3])/4])) *)  • That's right! In fact I've chosen this example because the exact solution can be found but forgot to write about it. Consider the same integral with the argument of the first Bessel function now being$K_0(|x_1-x_3|)$. The substitution trick becomes more problematic but perhaps doable. Now consider where I have say four variables$x_i$and the number of pairwise differences$x_i-x_j$(6) is larger then the number of variables (4). I do not see a way to bypass the numeric integration in this case. – Weather Report Oct 26 '20 at 19:26 • This is why it's important to ask about your actual question, not a substitute. – Roman Oct 26 '20 at 19:27 • I need to compute numerous integrals of this form and this particular one I thought captured the problem. Your point is taken and your effort appreciated. I actually might reconsider my computational approach based on your suggestion. – Weather Report Oct 26 '20 at 19:38 • @Roman What is the relation between the Integrate you give in the first codeline and the first codeline NIntegrate @Weather Report asks for? – Ulrich Neumann Oct 26 '20 at 21:53 • @Roman Thanks a lot. Your answer is very impressive for me, now trying to understand the substitution x1,x2,x3-> r1,r2,r3 in detail. – Ulrich Neumann Oct 27 '20 at 7:09 final modification ;-) Without "Subscripts" ( avoid use of subscripts ) realabs = Sqrt[#^2] & NIntegrate[r1 r2 r3 BesselK[0, realabs[r1]]^2 BesselK[0,Abs[r1 - E^(I \[Phi]2) r2]] BesselK[0,Abs[E^(I \[Phi]2) r2 - E^(I \[Phi]3) r3]]^3, {\[Phi]2, 0, 2 \[Pi]}, {\[Phi]3, 0, 2 \[Pi]}, {r1, 0, \[Infinity]}, {r2,0, \[Infinity]},{r3, 0, \[Infinity]} , WorkingPrecision -> 15 , AccuracyGoal -> 5,Method -> {"AdaptiveMonteCarlo","SymbolicProcessing" -> 0}] // AbsoluteTiming (*{60.5697, 7.37415114501064}*)  evalutes the intergal without error messages! • Unfortunately I do not have the latest Mathematica (v.11 here) and it seems that RealAbs is not available to me. However I do not understand (1) why do you suggest to use RealAbs on a complex-valued expressions? (2) How can the integral be zero given that the integrand is strictly positive? – Weather Report Oct 26 '20 at 16:21 • @WeatherReport RealAbs was introduced , because Abs sometimes shows problems with NDSolve, NIntegrate,.... The parameters r1,r2,r3  are real, that's why RealAbs is permitted. Try as a workaround realabs= Sqrt[#^2]&. – Ulrich Neumann Oct 26 '20 at 16:32 • Sorry I'm not following. I use genuinely complex coordinates to parametrize 2d plane so your suggestion should not work directly. I have tried to use realabs[z_] := Sqrt[Re[z Conjugate[z]]] without any change. It is very suspicious to me that you get zero (the expression in integral is positive!) and that exact integration can compute anything at all in this case. – Weather Report Oct 26 '20 at 16:55 • @WeatherReport See the modification in my answer, hopefully that's the solution you're looking for! – Ulrich Neumann Oct 26 '20 at 18:27 • Thank you, this seems to work better! However as far as I can tell the accuracy goal is far from reached. On three different runs I've got 7.07916951265200, 6.38668393270717 and 7.50130413248911. I actually have an impression that Monte-Carlo methods do not care too much about the accuracy goals:) – Weather Report Oct 26 '20 at 19:33 We can try to use compiled function, then we have f = Compile[{{r1, _Real}, {r2, _Real}, {r3, _Real}, {p2, _Real}, {p3, \ _Real}}, r1 r2 r3 BesselK[0, r1]^2 BesselK[0, Abs[r1 - E^(I p2) r2]] BesselK[0, Abs[E^(I p2) r2 - E^(I p3) r3]]^3, Parallelization -> True]; zero =$MachineEpsilon; L = Infinity;
NIntegrate[
f[r1, r2, r3, p2, p3], {r1, zero, L}, {r2, zero, L}, {r3, zero,
L}, {p2, 0, 2 \[Pi]}, {p3, 0, 2 \[Pi]},
Method -> "QuasiMonteCarlo"] // Quiet // AbsoluteTiming
{1.34761, 11.3053}


It is much quicker than other options and also close to the exact solution. Another method even faster but it fluctuates much

NIntegrate[
f[r1, r2, r3, p2, p3], {r1, zero, L}, {r2, zero, L}, {r3, zero,
L}, {p2, 0, 2 \[Pi]}, {p3, 0, 2 \[Pi]}, Method -> "MonteCarlo"] //
Quiet // AbsoluteTiming

(*Out[]= {1.24183, 11.4401}*)

• Hi, thanks for the suggestion! I've tested it on another example and unfortunately it does not perform as well. I'll update the question. – Weather Report Nov 8 '20 at 12:22
• @WeatherReport It could be better you show all integrals you try to calculate. There is no general algorithm for any case, but we can do something in every particular case. – Alex Trounev Nov 8 '20 at 12:39
• So your approach with the compiled functions turned out to be very helpful. Then it turns out I don't even have to use the MonteCarlo methods and standard GlobalAdaptive works fine and more reliably. Another thing that helped me is a realization that I do not have to consider all integrals one-by-one (I need only their sum) but instead compute only the total integral. So now I do not need to optimize every piece separately. All answers in this thread are very helpful, but yours nailed down my problem. Thanks! – Weather Report Nov 10 '20 at 10:12
• @WeatherReport You are welcome! – Alex Trounev Nov 10 '20 at 10:38