1
$\begingroup$

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?

$\endgroup$
1
  • 2
    $\begingroup$ 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 ... $\endgroup$ Jun 18 '20 at 19:34
3
$\begingroup$

A direct approach, bypassing NIntegrate:

(* Gauss-Kronrod rule abscissae, weights, and error weights *)
{abs, wts, err} = 
  Developer`ToPackedArray /@ 
   NIntegrate`GaussKronrodRuleData[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  *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.