How to calculate the double integral of a double integral?

Question

May I ask is there any wise way to calculate the following integration?

$\int_{-1}^{1}\int_{-1}^{1}\int_{-1}^{1}\int_{-1}^{1}\frac{\left( 0.25 \left(-1+\eta _1\right) \eta _1 \left(-1+2 \eta _2\right)\right) \left( \xi _1- \eta _1\right){}^2 \left(-0.25 \left(-1+\xi _1\right) \xi _1 \left(-1+2 \xi _2\right)\right)}{8 \left(4+\left(\xi _1- \eta _1\right){}^2+\left(\xi _2- \eta _2\right){}^2\right){}^{3/2}} d\eta_2 d\eta_1 d\xi_2 d\xi_1$

I have tried to use command $Int$ in Mathematica but it is not likely to yield the solution.

So far, I have tried to calculate the integral separately as follows

$\text{temp1} = \text{Integrate}\left[-\frac{\left(-0.2 \left(-1+\eta _1\right) \eta _1 \left(-1+2 \eta _2\right)\right) \left(-\eta _1+\xi _1\right){}^2 \left(-0.25 \left(-1+\xi _1\right) \xi _1 \left(-1+2 \xi _2\right)\right)}{8 \left(4+\left(-\eta _1+\xi _1\right){}^2+\left(- \eta _2+\xi _2\right){}^2\right){}^{3/2}},\left\{\eta _2,-1,1\right\},\left\{\eta _1,-1,1\right\}\right]$

$\text{temp2} = \text{Integrate}\left[\text{temp1},\left\{\xi _2,-1,1\right\},\left\{\xi _1,-1,1\right\}\right]$

I have also tried to integrate with $\eta_1$ only, and the simplify to get the result $-\frac{0.0024868 \left(\eta _1-1\right) \eta _1 \left(\xi _1-1\right) \xi _1 \left(2 \xi _2-1\right) \left(1. \eta _1-1. \xi _1\right){}^2 \left(\frac{4. \eta _1 \xi _1-2. \eta _1^2-2. \xi _1^2-2. \xi _2^2+3. \xi _2-9.}{\sqrt{-2. \eta _1 \xi _1+1. \eta _1^2+1. \xi _1^2+1. \xi _2^2-2. \xi _2+5.}}+\frac{-4. \eta _1 \xi _1+2. \eta _1^2+2. \xi _1^2+2. \xi _2^2+1. \xi _2+7.}{\sqrt{-2. \eta _1 \xi _1+1. \eta _1^2+1. \xi _1^2+1. \xi _2^2+2. \xi _2+5.}}\right)}{-2. \eta _1 \xi _1+1. \eta _1^2+1. \xi _1^2+4.}$

But after this, I cannot go any further.

Thank you.

Answer 1

This is not an answer but too long for a comment. The integral although may appear innocent is very nasty, I suggest, try Apart on it and you will see the source of the problems, if you try to do just one integral,

Integrate[((η1 - 1)*η1*(2*η2 - 
    1)*(ξ1 - η1)^2*(ξ1 - 1)*(2*ξ2 - 
    1))/(4 + (ξ1 - η1)^2 + (ξ2 - η2)^2)^(3/2), {\
η1, -1, 1}, 
 Assumptions -> {-1 < η2 < 1, -1 < ξ1 < 1, -1 < ξ2 < 1}]

and have enough patience, you will get a few lines' result, which you need to integrate again, get something even nastier and so on. I see two ways around it, a clever manipulation and substitution, there may be a way with some trigonometric substitution but this will require probably a lot of time if you don't do it on a regular basis.

The other way that I can think of is NIntegrate, I tried that but Mathematica is unhappy with that reporting slow convergence but eventually spits out a results.

A quick edit, what I get after integrating just the code above is this:

1/2 (-1 + 2 η2) (-1 + ξ1) (-1 + 
   2 ξ2) ((-5 Sqrt[
       4 + η2^2 + ξ1^2 - 
        2 η2 ξ2 + ξ2^2] - η2^2 Sqrt[
       4 + η2^2 + ξ1^2 - 2 η2 ξ2 + ξ2^2] + 
      27 ξ1 Sqrt[
       4 + η2^2 + ξ1^2 - 2 η2 ξ2 + ξ2^2] + 
      5 η2^2 ξ1 Sqrt[
       4 + η2^2 + ξ1^2 - 2 η2 ξ2 + ξ2^2] - 
      11 ξ1^2 Sqrt[
       4 + η2^2 + ξ1^2 - 2 η2 ξ2 + ξ2^2] + 
      5 ξ1^3 Sqrt[
       4 + η2^2 + ξ1^2 - 2 η2 ξ2 + ξ2^2] + 
      2 η2 ξ2 Sqrt[
       4 + η2^2 + ξ1^2 - 2 η2 ξ2 + ξ2^2] - 
      10 η2 ξ1 ξ2 Sqrt[
       4 + η2^2 + ξ1^2 - 
        2 η2 ξ2 + ξ2^2] - ξ2^2 Sqrt[
       4 + η2^2 + ξ1^2 - 2 η2 ξ2 + ξ2^2] + 
      5 ξ1 ξ2^2 Sqrt[
       4 + η2^2 + ξ1^2 - 2 η2 ξ2 + ξ2^2] + 
      16 Sqrt[5 + η2^2 - 2 ξ1 + ξ1^2 - 
        2 η2 ξ2 + ξ2^2] + 
      4 η2^2 Sqrt[
       5 + η2^2 - 2 ξ1 + ξ1^2 - 
        2 η2 ξ2 + ξ2^2] - 
      20 ξ1 Sqrt[
       5 + η2^2 - 2 ξ1 + ξ1^2 - 
        2 η2 ξ2 + ξ2^2] - 
      5 η2^2 ξ1 Sqrt[
       5 + η2^2 - 2 ξ1 + ξ1^2 - 
        2 η2 ξ2 + ξ2^2] + 
      4 ξ1^2 Sqrt[
       5 + η2^2 - 2 ξ1 + ξ1^2 - 
        2 η2 ξ2 + ξ2^2] - 
      5 ξ1^3 Sqrt[
       5 + η2^2 - 2 ξ1 + ξ1^2 - 
        2 η2 ξ2 + ξ2^2] - 
      8 η2 ξ2 Sqrt[
       5 + η2^2 - 2 ξ1 + ξ1^2 - 
        2 η2 ξ2 + ξ2^2] + 
      10 η2 ξ1 ξ2 Sqrt[
       5 + η2^2 - 2 ξ1 + ξ1^2 - 
        2 η2 ξ2 + ξ2^2] + 
      4 ξ2^2 Sqrt[
       5 + η2^2 - 2 ξ1 + ξ1^2 - 
        2 η2 ξ2 + ξ2^2] - 
      5 ξ1 ξ2^2 Sqrt[
       5 + η2^2 - 2 ξ1 + ξ1^2 - 
        2 η2 ξ2 + ξ2^2] + 
      Sqrt[(4 + η2^2 + ξ1^2 - 
          2 η2 ξ2 + ξ2^2) (5 + η2^2 - 
          2 ξ1 + ξ1^2 - 
          2 η2 ξ2 + ξ2^2)] (3 η2^2 + 2 ξ1 - 
         2 ξ1^2 - 6 η2 ξ2 + 
         3 (4 + ξ2^2)) Log[-ξ1 + Sqrt[
         4 + η2^2 + ξ1^2 - 2 η2 ξ2 + ξ2^2]] - 
      Sqrt[(4 + η2^2 + ξ1^2 - 
          2 η2 ξ2 + ξ2^2) (5 + η2^2 - 
          2 ξ1 + ξ1^2 - 
          2 η2 ξ2 + ξ2^2)] (3 η2^2 + 2 ξ1 - 
         2 ξ1^2 - 6 η2 ξ2 + 3 (4 + ξ2^2)) Log[
        1 - ξ1 + Sqrt[
         5 + η2^2 - 2 ξ1 + ξ1^2 - 
          2 η2 ξ2 + ξ2^2]])/(Sqrt[(4 + η2^2 + \
ξ1^2 - 2 η2 ξ2 + ξ2^2) (5 + η2^2 - 
        2 ξ1 + ξ1^2 - 2 η2 ξ2 + ξ2^2)]) + (31 + 
      7 η2^2 - 15 ξ1 - 5 η2^2 ξ1 - 3 ξ1^2 - 
      5 ξ1^3 - 14 η2 ξ2 + 10 η2 ξ1 ξ2 + 
      7 ξ2^2 - 5 ξ1 ξ2^2 - 
      4 Sqrt[(4 + η2^2 + ξ1^2 - 
          2 η2 ξ2 + ξ2^2) (5 + η2^2 + 
          2 ξ1 + ξ1^2 - 2 η2 ξ2 + ξ2^2)] + 
      5 ξ1 Sqrt[(4 + η2^2 + ξ1^2 - 
          2 η2 ξ2 + ξ2^2) (5 + η2^2 + 
          2 ξ1 + ξ1^2 - 2 η2 ξ2 + ξ2^2)] + 
      Sqrt[5 + η2^2 + 2 ξ1 + ξ1^2 - 
        2 η2 ξ2 + ξ2^2] (3 η2^2 + 2 ξ1 - 
         2 ξ1^2 - 6 η2 ξ2 + 
         3 (4 + ξ2^2)) Log[ξ1 + Sqrt[
         4 + η2^2 + ξ1^2 - 2 η2 ξ2 + ξ2^2]] - 
      Sqrt[5 + η2^2 + 2 ξ1 + ξ1^2 - 
        2 η2 ξ2 + ξ2^2] (3 η2^2 + 2 ξ1 - 
         2 ξ1^2 - 6 η2 ξ2 + 3 (4 + ξ2^2)) Log[
        1 + ξ1 + Sqrt[
         5 + η2^2 + 2 ξ1 + ξ1^2 - 
          2 η2 ξ2 + ξ2^2]])/(Sqrt[
     5 + η2^2 + 2 ξ1 + ξ1^2 - 
      2 η2 ξ2 + ξ2^2]))

It is interesting that you get such a simple result by that integration.

Answer 2

If "is there any wise way?" doesn't exclude a Monte Carlo method, then the following might be considered:

f = ((1/4) (n1 - 1) n1 (2 n2 - 1)) (e1 - n1)^2 (-(1/4) (e1 - 1) e1 (2 e2 - 1))/
    (8 (4 + (e1 - n1)^2 + (e2 - n2)^2)^(3/2));
NIntegrate[f, {e1, -1, 1}, {e2, -1, 1}, {n1, -1, 1}, {n2, -1, 1}, 
    Method -> "AdaptiveQuasiMonteCarlo"]
(* 0.0007126203438669235 *)

Update

Someone more experienced with tweaking the options of NIntegrate is needed. I tried the following brute force method and obtained a different enough answer to be concerned about just using the default settings for NIntegrate:

(* Generate a bunch of random numbers between -1 and +1 *)
n = 100000000;
e1 = RandomReal[{-1, 1}, n];
e2 = RandomReal[{-1, 1}, n];
n1 = RandomReal[{-1, 1}, n];
n2 = RandomReal[{-1, 1}, n];

(* Estimate of integral *)
16 Mean[((1/4) (n1 - 1) n1 (2 n2 - 1)) (e1 - n1)^2 (-(1/4) (e1 - 1) e1 (2 e2 - 1))/
(8 (4 + (e1 - n1)^2 + (e2 - n2)^2)^(3/2))]
(* 0.0007193399714451682 *)

(* Standard error *)
16 StandardDeviation[((1/4) (n1 - 1) n1 (2 n2 - 1)) (e1 - n1)^2 (-(1/4) (e1 - 1) e1 (2 e2 - 1))/
(8 (4 + (e1 - n1)^2 + (e2 - n2)^2)^(3/2))]/n
(* 4.7445218486769696`*^-11 *)

Adding some tweaks to NIntegrate I get a result very close to the brute force method but no meeting the convergence criteria:

f = ((1/4) (n1 - 1) n1 (2 n2 - 1)) (e1 - n1)^2 (-(1/4) (e1 - 1) e1 (2 e2 - 1))/
  (8 (4 + (e1 - n1)^2 + (e2 - n2)^2)^(3/2));
NIntegrate[f, {e1, -1, 1}, {e2, -1, 1}, {n1, -1, 1}, {n2, -1, 1}, 
  Method -> "AdaptiveQuasiMonteCarlo", MinRecursion -> 9, 
  MaxPoints -> 100, WorkingPrecision -> 10]
(* 0.00071934046972184822351421014894205721`10. *)

NIntegrate::maxp: The integral failed to converge after 25600 integrand evaluations. NIntegrate obtained 0.0007193404697 and 0.00001781325045. for the integral and error estimates. >>

Answer 3

Let us use 18 methods of numerical integration and see that the results of which overlaps.

Warning!!!.The calculations take a long time. (several hours).

  ClearAll["Global`*"]

  f = ((1/4) (n1 - 1) n1 (2 n2 - 1)) (e1 - n1)^2 (-(1/4) (e1 - 1) e1 
  (2 e2 - 1))/(8 (4 + (e1 - n1)^2 + (e2 - n2)^2)^(3/2));

  method = {"GlobalAdaptive", "LocalAdaptive", 
  "Trapezoidal", "MultiPeriodic", "MonteCarlo", "QuasiMonteCarlo", 
  "AdaptiveMonteCarlo", "DoubleExponentialOscillatory", 
  "TrapezoidalRule", "NewtonCotesRule", 
  "GaussBerntsenEspelidRule", "GaussKronrodRule", 
  "LobattoKronrodRule", "ClenshawCurtisRule", "MultipanelRule", 
  "CartesianRule", "MultidimensionalRule", "MonteCarloRule"};

  Grid[Prepend[Transpose[{method, NIntegrate[f, {e1, -1, 1}, {e2, -1, 1}, {n1, -1, 1}, {n2, -1, 1}, 
  Method -> #,
  Exclusions -> ((8 (4 + (e1 - n1)^2 + (e2 - n2)^2)^(3/2)) == 0)] & /@ method // Quiet}], 
  {"Method Name", "Result"}], 
  Background -> {None, {Lighter[Yellow, .9], {White, 
  Lighter[Blend[{Blue, Green}], .8]}}}, 
  Dividers -> {{Darker[Gray, .6], {Lighter[Gray, .5]}, 
  Darker[Gray, .6]}, {Darker[Gray, .6], Darker[Gray, .6], {False}, 
  Darker[Gray, .6]}}, Alignment -> {{Left, Left}}, 
  ItemSize -> {{20, 10}}, Frame -> Darker[Gray, .6], ItemStyle -> 14, 
  Spacings -> {Automatic, .8}, Frame -> All]

The most likely result is 0.000719295