Nested integral takes too long with NIntegrate

I am attempting to calculate a nested multiple integral. The code snippet is given by

 i1[x1_, x2_, x3_, x4_?NumericQ] := i1[x1, x2, x3, x4] =
NIntegrate[(x1 x2 x3 x4 x5 (1 - x1 - x2 - x3 - x4 - x5))^4/(
x2 x3 x4 x5 (1 - x1 - x2 - x3 - x4 - x5) +
x1 x3 x4 x5 (1 - x1 - x2 - x3 - x4 - x5) +
x1 x2 x4 x5 (1 - x1 - x2 - x3 - x4 - x5) +
x1 x2 x3 x5 (1 - x1 - x2 - x3 - x4 - x5) +
x1 x2 x3 x4 (1 - x1 - x2 - x3 - x4 - x5) +
x1 x2 x3 x4 x5)^6, {x5, 0, 1 - x1 - x2 - x3 - x4}];

i2[x1_, x2_, x3_?NumericQ] :=i2[x1, x2, x3] =
NIntegrate[i1[x1, x2, x3, x4], {x4, 0, 1 - x1 - x2 - x3}];

i3[x1_, x2_?NumericQ] := i3[x1, x2] = NIntegrate[i2[x1, x2, x3], {x3, 0, 1 - x1 - x2}];

i4[x1_?NumericQ] := i4[x1] = NIntegrate[i3[x1, x2], {x2, 0, 1 - x1}];
NIntegrate[i4[x1], {x1, 0, 1}]


However this seems to take way too long (I even left it overnight and it had not been evaluated). Can anyone help me to get around this issue?

• This NIntegrate[(x1 x2 x3 x4 x5 (1-x1-x2-x3-x4-x5))^4/(x2 x3 x4 x5 (1-x1-x2-x3-x4-x5)+x1 x3 x4 x5 (1-x1-x2-x3-x4-x5)+x1 x2 x4 x5 (1-x1-x2-x3-x4-x5)+x1 x2 x3 x5 (1-x1-x2-x3-x4-x5)+x1 x2 x3 x4 (1-x1-x2-x3-x4-x5)+x1 x2 x3 x4 x5)^6,{x1,0,1},{x2,0,1-x1},{x3,0,1-x1-x2},{x4,0,1-x1-x2-x3},{x5,0,1-x1-x2-x3-x4},Method->"MultidimensionalRule",PrecisionGoal->3] gives an answer, but also error messages ... Jun 18 '20 at 19:34

A direct approach, bypassing NIntegrate:

(* Gauss-Kronrod rule abscissae, weights, and error weights *)
{abs, wts, err} =
DeveloperToPackedArray /@
NIntegrateGaussKronrodRuleData[21, MachinePrecision];

(* integrand evaluations *)
data =
With[{i = (
(x1 x2 x3 x4 x5 (1 - x1 - x2 - x3 - x4 - x5))^4/
(x2 x3 x4 x5 (1 - x1 - x2 - x3 - x4 - x5) +
x1 x3 x4 x5 (1 - x1 - x2 - x3 - x4 - x5) +
x1 x2 x4 x5 (1 - x1 - x2 - x3 - x4 - x5) +
x1 x2 x3 x5 (1 - x1 - x2 - x3 - x4 - x5) +
x1 x2 x3 x4 (1 - x1 - x2 - x3 - x4 - x5) +
x1 x2 x3 x4 x5)^6 //
# (1 - x1) (1 - x1 - x2) (1 - x1 - x2 - x3) *
(1 - x1 - x2 - x3 - x4) &
) /. x5 -> xj (1 - x1 - x2 - x3 - x4) /.
x4 -> x4 (1 - x1 - x2 - x3) /.
x3 -> x3 (1 - x1 - x2) /.
x2 -> x2 (1 - x1)},
Compile[{{xj, _Real, 1}, {x4, _Real}},
Table[i, {x1, xj}, {x2, xj}, {x3, xj}],
RuntimeAttributes -> {Listable}, Parallelization -> True
][abs, abs]
]; // AbsoluteTiming
(*  {4.68719, Null}  *)

(* integral *)
Nest[#.wts &, data, 5]
(*  0.00835538  *)

(* error estimate: ~4 digits of precision *)
Nest[#.err &, data, 5]
(*  3.43986*10^-6  *)