# How to calculate the double integral of a double integral?

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.

• Integrate[f,{x,xmin,xmax},{y,ymin,ymax},.......] gives the multiple integral. Post the code you've tried so that we can see what is wrong. There's no limit on the number of successive integrals you can define – Jason B. Apr 18 '16 at 12:30
• Note that Mathematica treats 0.25 as a machine-precision real number, whereas it would treat 1/4 as an exact rational number. See here for further information on the distinction. Mathematica sometimes uses different algorithms to treat rational vs. real numbers, so may have some success if you use only rational (or only real) numbers. – Michael Seifert May 18 '16 at 16:23

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.

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},
(* 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.0007193404697218482235142101489420572110. *)


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

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));

"Trapezoidal", "MultiPeriodic", "MonteCarlo", "QuasiMonteCarlo",
"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

• I like it. Is the term "most likely result" a statistical or engineering term? – JimB May 18 '16 at 17:58
• @JimBaldwin.On faith term. – Mariusz Iwaniuk May 18 '16 at 19:58