# Help in solving high-dimensional definite integral

I am trying to solve a 6D integral over this integrated

integrand=ConditionalExpression[(9*(1 + E^(R*\[Eta]1))^2*(1 +E^(R*\[Eta]2))^2*R^6*(\[Eta]1 - \[Xi]1)*(\[Eta]1 + \[Xi]1)*(\[Eta]2 - \[Xi]2)*(\[Eta]2 + \[Xi]2))/(E^(R*(-2 + \[Eta]1 + \[Eta]2 + \[Xi]1 + \[Xi]2))*(128*Pi^2*(3 + 3*E^R + R*(3 + R))^2*Sqrt[R^2*(-2 + \[Eta]1^2 + \[Eta]2^2 + \[Xi]1^2 - 2*\[Eta]1*\[Eta]2*\[Xi]1*\[Xi]2 + \[Xi]2^2 - 2*Sqrt[-((-1 + \[Eta]1^2)*(-1 + \[Xi]1^2))]*Sqrt[-((-1 + \[Eta]2^2)*(-1 + \[Xi]2^2))]*Cos[\[Phi]1 - \[Phi]2])])), Re[R] > 0]

Vee = Integrate[integrand, {\[Xi]1, 1, \[Infinity]}, {\[Eta]1, -1, 1}, {\[Phi]1, 0, 2*\[Pi]}, {\[Xi]2, 1, \[Infinity]}, {\[Eta]2, -1, 1}, {\[Phi]2, 0, 2*\[Pi]}]


with the limits $$\xi_1,\xi_2$$ from 1 to $$\infty$$, $$\eta_1$$,$$\eta_2$$ from -1 to 1 and $$\phi_1,\phi_2$$ from 0 to $$2\pi$$. By the way, this is calculating the electrostatic energy of the H$$_2$$ molecule in confocal, elliptic coordinates. The problem is, this integral is running for a long time and does not finish. I am not sure if it is working properly or not, since there is no response, it just keeps running. Any idea how to solve or speed up this integral?

• People here generally like users to post code as copyable Mathematica code instead of just images or TeX, so they can copy-paste it. It makes it convenient for them and more likely you will get someone to help you. You may find the meta Q&A, How to copy code from Mathematica so it looks good on this site, helpful May 22, 2023 at 0:47
• I tried to follow the procedure. Is it better now? May 22, 2023 at 0:54
• yes, thanks. Not sure whether it will be easy for someone to help or not. If the conditional expression is meant to add an assumption that Re[R] > 0, the typical way to do that is adding the option Assumptions -> Re[R] > 0 to Integrate. May 22, 2023 at 1:01
• I will try that, definitely R is a real number and larger than 0 May 22, 2023 at 1:02
• "R is a real number and larger than 0" -- then use Assumptions -> R > 0. Variables that appear in inequalities are assumed to be real. Re[R] > 0 signals that R is complex, with a positive real part. (Integrate[] still seems to run forever, though.) May 22, 2023 at 1:05

We can reduce dimension from 6 to 4 using

 Integrate[
1/Sqrt[a Cos[p1 - p2] + b], {p1, 0, 2 Pi}, {p2, 0, 2 Pi},
Assumptions -> {a > 0, b > 0}]

(*Out[1]= (8 \[Pi] EllipticK[(2 a)/(a + b)])/Sqrt[a + b] *)
% /. {a -> -2*Sqrt[-((\[Eta]1^2 - 1)*(\[Xi]1^2 - 1))]*Sqrt[-((\[Eta]2^2 - 1)*(\[Xi]2^2 - 1))], b -> \[Eta]1^2 - 2*\[Eta]1*\[Eta]2*\[Xi]1*\[Xi]2 + \[Eta]2^2 + \[Xi]1^2 + \[Xi]2^2 - 2}

(*Out[2]= (8 \[Pi] EllipticK[-((
4 Sqrt[-((-1 + \[Eta]1^2) (-1 + \[Xi]1^2))]
Sqrt[-((-1 + \[Eta]2^2) (-1 + \[Xi]2^2))])/(-2 + \[Eta]1^2 + \
\[Eta]2^2 + \[Xi]1^2 - 2 \[Eta]1 \[Eta]2 \[Xi]1 \[Xi]2 + \[Xi]2^2 -
2 Sqrt[-((-1 + \[Eta]1^2) (-1 + \[Xi]1^2))]
Sqrt[-((-1 + \[Eta]2^2) (-1 + \[Xi]2^2))]))])/Sqrt[-2 + \
\[Eta]1^2 + \[Eta]2^2 + \[Xi]1^2 -
2 \[Eta]1 \[Eta]2 \[Xi]1 \[Xi]2 + \[Xi]2^2 -
2 Sqrt[-((-1 + \[Eta]1^2) (-1 + \[Xi]1^2))]
Sqrt[-((-1 + \[Eta]2^2) (-1 + \[Xi]2^2))]]*)


Therefore integral can be evaluated in 4 dimensional space as follows

f[R_?NumericQ] := (9*R^5)/(128*
Pi^2*(R*(R + 3) + 3*E^R +
3)^2) NIntegrate[(8*(1 + E^(R*\[Eta]1))^2*(1 +
E^(R*\[Eta]2))^2*
Pi*(\[Eta]1 - \[Xi]1)*(\[Eta]1 + \[Xi]1)*(\[Eta]2 - \[Xi]2)*(\
\[Eta]2 + \[Xi]2)*
EllipticK[-((4*Sqrt[-((-1 + \[Eta]1^2)*(-1 + \[Xi]1^2))]*
Sqrt[-((-1 + \[Eta]2^2)*(-1 + \[Xi]2^2))])/(-2 + \
\[Eta]1^2 + \[Eta]2^2 + \[Xi]1^2 -
2*\[Eta]1*\[Eta]2*\[Xi]1*\[Xi]2 + \[Xi]2^2 -
2*Sqrt[-((-1 + \[Eta]1^2)*(-1 + \[Xi]1^2))]*
Sqrt[-((-1 + \[Eta]2^2)*(-1 + \[Xi]2^2))]))])/(E^(R*(-2 \
+ \[Eta]1 + \[Eta]2 + \[Xi]1 + \[Xi]2))*
Sqrt[-2 + \[Eta]1^2 + \[Eta]2^2 + \[Xi]1^2 -
2*\[Eta]1*\[Eta]2*\[Xi]1*\[Xi]2 + \[Xi]2^2 -
2*Sqrt[-((-1 + \[Eta]1^2)*(-1 + \[Xi]1^2))]*
Sqrt[-((-1 + \[Eta]2^2)*(-1 + \[Xi]2^2))]]), {\[Xi]1, 1,
Infinity}, {\[Eta]1, -1, 1}, {\[Xi]2, 1, Infinity}, {\[Eta]2, -1,


Example of usage

f[1] // AbsoluteTiming

Out[]= {0.412131, 0.588934}


Also this integral can be computed in 6D as follows

integrand = (9*R^6*(\[Eta]1 - \[Xi]1)*(\[Eta]1 + \[Xi]1)*(\[Eta]2 - \[Xi]2)*(\[Eta]2 + \[Xi]2)*(E^(\[Eta]1*R) + 1)^2*
(E^(\[Eta]2*R) + 1)^2)/(E^(R*(\[Eta]1 + \[Eta]2 + \[Xi]1 + \[Xi]2 - 2))*(128*Pi^2*(R*(R + 3) + 3*E^R + 3)^2*
Sqrt[R^2*(-(2*Sqrt[-((\[Eta]1^2 - 1)*(\[Xi]1^2 - 1))]*Sqrt[-((\[Eta]2^2 - 1)*(\[Xi]2^2 - 1))]*
Cos[\[Phi]1 - \[Phi]2]) + \[Eta]1^2 - 2*\[Eta]1*\[Eta]2*\[Xi]1*\[Xi]2 + \[Eta]2^2 + \[Xi]1^2 + \[Xi]2^2 - 2)]));
lst = Table[{x,
NIntegrate[
integrand /. R -> x, {\[Xi]1, 1, Infinity}, {\[Eta]1, -1,
1}, {\[Phi]1, 0, 2*Pi}, {\[Xi]2, 1, Infinity}, {\[Eta]2, -1,
1}, {\[Phi]2, 0, 2*Pi},
Method -> "AdaptiveQuasiMonteCarlo"]}, {x, 1/10, 5, 1/10}]


Visualization

ListPlot[{lst, {{1, f[1]}}}, PlotStyle -> {Blue, Red}, Frame -> True,
FrameLabel -> {"R", "Vee"}, PlotLegends -> {"6D", "4D"}]


• Yeah I guess numerical integration is better in this case. I was hoping for an analytic solution, but Mathematica seems not to be able to give it to me. May 31, 2023 at 0:34
• @Guiste Yes, you are right. Analytical expression in this case maybe is not so useful. Even reduced from 6D in 4D looks very cumbersome. May 31, 2023 at 4:25