How to make this kind of integration with the trigonometric function?

Question

I am trying to make the following integration. But I am not sure about the result obtained as follows.

Is there some way to get more reliable and accurate results?

Thanks.

NIntegrate[
    (x1 x2 x3)/((1 + x1^2)^(3/2) (1 + x3^2)^(3/2) (1 + x2^2 - 7 x2 Cos[x4]) Sqrt[x1^2 + x2^2 - x1 x2 Cos[x5]] Sqrt[ 1 + 9 (x1^2 + x2^2 - x1 x2 Cos[x6])])
     , {x1, 0, ∞}
     , {x2, 0, ∞}
     , {x3, 0, ∞}
     , {x4, 0, 2 π}
     , {x5, 0, 2 π}
     , {x6, 0, 2 π}
]

Edit:

To avoid the singularity problem, I revise the above problem and make the integral as follows,

        fexp = (x1 x2 x3)/((1 + x1^2)^(3/2) (1 + x3^2)^(
        3/2) (1 + x2^2 - x2 Cos[x4]) Sqrt[
        3 + x1^2 + x2^2 - 2 x1 x2 Cos[x5]] Sqrt[
        1 + x1^2 + x2^2/4 - x1 x2 Cos[x6]]);

    NIntegrate[fexp
               , {x1, 0, ∞}
               , {x2, 0, ∞}
               , {x3, 0, ∞}
               , {x4, 0, 2 π}
               , {x5, 0, 2 π}
               , {x6, 0, 2 π}]

   Out=75.0157  (with warning messages, Numerical integration converging 
                too slowly; .... )

In addition, I want to get the results numerically by taking some proper algorithms. Of course, we can check the results analytically (partly).

Answer 1

You can do most of the inegrations analyticlly. As the other comments mention, you have to consider singularities.

f = (x1 x2 x3)/((1 + x1^2)^(3/2) (1 + x3^2)^(3/2) (1 + x2^2 - 
        7 x2 Cos[x4]) Sqrt[x1^2 + x2^2 - x1 x2 Cos[x5]] Sqrt[
       1 + 9 (x1^2 + x2^2 - x1 x2 Cos[x6])]) // Together;

int21 = Integrate[f, {x5, 0, 2 \[Pi]}, {x6, 0, 2 \[Pi]}, 
  Assumptions -> {x1 > 0, x2 > 0, x3 > 0, 0 < x4 < 2 Pi}]

(*   (16 x1 x2 x3 EllipticK[-((2 x1 x2)/(
     x1^2 - x1 x2 + x2^2))] EllipticK[-((18 x1 x2)/(
     1 + 9 (x1^2 - x1 x2 + x2^2)))])/(Sqrt[(x1^2 - x1 x2 + x2^2) (1 + 
      9 (x1^2 - x1 x2 + x2^2))] ((1 + x1^2) (1 + x3^2))^(
   3/2) (1 + x2^2 - 7 x2 Cos[x4]))   *)

int22 = Integrate[int21, {x3, 0, Infinity}, 
  Assumptions -> {x1 > 0, x2 > 0, x3 > 0, 0 < x4 < 2 Pi}]

(*   (16 x1 x2 EllipticK[-((2 x1 x2)/(x1^2 - x1 x2 + x2^2))] EllipticK[-((
   18 x1 x2)/(1 + 9 (x1^2 - x1 x2 + x2^2)))])/((1 + x1^2)^(
 3/2) Sqrt[(x1^2 - x1 x2 + x2^2) (1 + 9 (x1^2 - x1 x2 + x2^2))] (1 + 
   x2^2 - 7 x2 Cos[x4]))   *)

Singularities at

ContourPlot[(1 + x2^2 - 7 x2 Cos[x4]) == 0, {x2, 0, 10}, {x4, 0, 
  2 Pi}]

In this range, integral is zero.

int23 = Integrate[int22, {x4, 0, 2 Pi}, 
  Assumptions -> x1 > 0 && x2 > 0, PrincipalValue -> True]

(*   ConditionalExpression[0, 
 3 Sqrt[5] + 2 x2 >= 7 && 2 x2 <= 7 + 3 Sqrt[5]]   *)

{red = Reduce[int23[[2]], x2], red // N}

(*   {1/2 (7 - 3 Sqrt[5]) <= x2 <= 1/2 (7 + 3 Sqrt[5]), 
 0.145898 <= x2 <= 6.8541}   *)

int24 = Integrate[int22, {x4, 0, 2 Pi}, 
 Assumptions -> 
  x1 > 0 && (0 < x2 < 1/2 (7 - 3 Sqrt[5]) || 
     1/2 (7 + 3 Sqrt[5]) < x2)]

(*   (32 \[Pi] x1 x2 EllipticK[-((2 x1 x2)/(
   x1^2 - x1 x2 + x2^2))] EllipticK[-((18 x1 x2)/(
   1 + 9 (x1^2 - x1 x2 + x2^2)))])/((1 + x1^2)^(
 3/2) Sqrt[(x1^2 - x1 x2 + x2^2) (1 - 47 x2^2 + x2^4) (1 + 
    9 (x1^2 - x1 x2 + x2^2))])   *)

int25 = NIntegrate[
   int24, {x1, 0, Infinity}, {x2, 0, 1/2 (7 - 3 Sqrt[5])}, 
   WorkingPrecision -> 15] + 
  NIntegrate[
   int24, {x1, 0, Infinity}, {x2, 1/2 (7 + 3 Sqrt[5]), Infinity}, 
   WorkingPrecision -> 15]

(*   3.820016280311350   *)

Edit

What i did with int23, is regarding x1 and x2 as parameters and calculating the principal value for x4. This is affirmed by inserting parameters.

ii22[x1_, x2_, x4_] = int22; Manipulate[
 Plot[ii22[x1, x2, x4], {x4, 0, 2 Pi}, PlotPoints -> 50, 
  PlotRange -> 5], {{x1, 2}, 0, 10, 
  Appearance -> "Labeled"}, {{x2, 1/E}, 0, 10, 
  Appearance -> "Labeled"}]

int23b = Integrate[int22 /. {x2 -> 1/E}, {x4, 0, 2 Pi}, 
  PrincipalValue -> True, Assumptions -> x1 > 0]

(*   0   *)

Reduce[0 < x4 < 2 Pi && 0 == (1 + x2^2 - 7 x2 Cos[x4]) /. 
  x2 -> 1/E, x4, Reals]

(*   x4 == 2 \[Pi] - ArcCos[(1 + E^2)/(7 E)] || 
 x4 == ArcCos[(1 + E^2)/(7 E)]   *)

NIntegrate[
 int22 /. {x1 -> 2, x2 -> 1/E}, {x4, 0, ArcCos[(1 + E^2)/(7 E)], 
  2 Pi - ArcCos[(1 + E^2)/(7 E)], 2 Pi}, Method -> "PrincipalValue"]

(*   3.71023*10^-14   *)

I'm no pure mathematician to be able to proof this kind of proceeding. But otherwise you must state, integral is not integrable.