How to define such a complex region and NIntegrate over it?

Question

This may be in my urgent, I'm currently dealing with an integration problem over the support of high dimensional [0, 1] cubic. The integration function itself is not that complicated, you can assume that

F[p1_, p2_, p3_, p4_, p5_, p6_, p7_, p8_, p9_, p10_] := p1 + 2*p2 + 3*p3 + 4*p4 + 5*p5 + 6*p6 + 7*p7 + 8*p8 + 9*p9 + 10*p10

But the region over cubic for integration is defined by the union of the following two region, i.e. $R_1 \cup R_2$:

$R_1: \prod_{i=1}^{s} p_{(i)} \leq c_1$

$R_2: \prod_{i=1}^{s} \frac{p_{(i)} - p_{(s+1)}}{1 - p_{(s+1)}} \leq c_2$

where $p_{(1)} \geq p_{(2)} \geq ...... \geq p_{(10)}$ is the ordered value of $(p_1, p_2, ... , p_{10})$, $0 < p_i < 1$ and $2 \leq s \leq 8$. $c_1, c_2$ are rather small value, for example, $c_1 = 0.0025, c_2 = 10^{-6}$.

I'm a primary learner in Mathematica, and could define some simple region only (like sphere, rectangle), but the region $R_1 \cup R_2$ is more complex than I expected, so I'm urgent to know how to define the region that ruled by relationship of ordering values.

Thanks very much!!

Answer 1

This is a hard problem. There are several difficulties to overсome: to write down $p_{(1)} \geq p_{(2)} \geq ...... \geq p_{(10)}$ in the author's notation, to integrate numerically over a complicated set in 10 dimensions. First, the ten-dimensional hypercube can be split into 10! ten-dimensional simplices (e.g. see that thread). For example, in three dimensions a similar partition looks as

Show[RegionPlot3D[ 0 <= x && x <= y && y <= z && z <= 1, {x, -1/2, 3/2}, {y, -1/2,  3/2}, 
{z, -1/2, 3/2}, PlotPoints -> 100,  PlotStyle -> Directive[Red], Mesh -> None], 
RegionPlot3D[ 0 <= z && z <= y && y <= x && x <= 1, {x, -1/2, 3/2}, {y, -1/2,  3/2}, 
{z, -1/2, 3/2}, PlotPoints -> 100, PlotStyle -> Directive[Orange], Mesh -> None], 
RegionPlot3D[ 0 <= y && y <= x && x <= z && z <= 1, {x, -1/2, 3/2}, {y, -1/2, 3/2}, 
{z, -1/2, 3/2}, PlotPoints -> 100, PlotStyle -> Directive[Yellow], Mesh -> None], 
RegionPlot3D[ 0 <= x && x <= z && z <= y && y <= 1, {x, -1/2, 3/2}, {y, -1/2,   3/2}, 
{z, -1/2, 3/2}, PlotPoints -> 100, PlotStyle -> Directive[Green], Mesh -> None], 
RegionPlot3D[ 0 <= z && z <= x && x <= y && y <= 1, {x, -1/2, 3/2}, {y, -1/2,   3/2}, 
{z, -1/2, 3/2}, PlotPoints -> 100, PlotStyle -> Directive[Blue], Mesh -> None], 
 RegionPlot3D[ 0 <= y && y <= z && z <= x && x <= 1, {x, -1/2, 3/2}, {y, -1/2,  3/2},   
{z, -1/2, 3/2}, PlotPoints -> 100,PlotStyle -> Directive[Gray], Mesh -> None]]

and in two dimensions a similar partition consists of two triangles.

Second, the relations 0 <= p10 && p10 <= p9 && p9 <= p8 && p8 <= p7 && p7 <= p6 && p6 <= p5 && p5 <= p4 && p4 <= p3 && p3 <= p2 && p2 <= p1 && p1 <= 1 are valid on one of the simplices. In view of it the integral under consideration over this simplex is

Integrate[(p1+2*p2+3*p3+4*p4+5*p5+6*p6+7*p7+8*p8+9*p9+10*p10)*
Piecewise[{{1,  p1*p2 <= 1/400 || ((-p3 + p2)* (-p2 + p1))/((1 - p3)* (1 - p2)) <= 1/1000000}, {0, True}}],
 {p1, 0, 1}, {p2, 0, p1}, {p3, 0, p2},{p4, 0, p3}, {p5, 0, p4}, {p6, 0, p5}, {p7, 0, p6}, {p8, 0, p7}, {p9, 0, p8},{p10,0,p9}]

I prefer a piecewise continuous integrand than integration over a complicated set. Now let us consider another simplex of the partition. There $p_{(1)} \geq p_{(2)} \geq ...... \geq p_{(10)}$ is a permutation of $p1,p2,\dots,p10$ (see three dimensions as a model) and we come to the integral

Integrate[($p_{(1)}$+ 2*$p_{(2)}$ + 3*$p_{(3)}$ + 4*$p_{(4)}$ + 5*$p_{(5)}$ + 6*$p_{(6)}$ + 7*$p_{(7)}$ + 8*$p_{(8)}$ + 9*$p_{(9)}$ + 10*$p_{(10)}$)* Piecewise[{{1, $p_{(1)}$$p_{(2)}$ <= 1/400 || ((- $p_{(3)}$+$p_{(2)}$) (-$p_{(2)}$ + $p_{(1)}$))/((1 - $p_{(3)}$)* (1 - $p_{(2)}$)) <= 1/1000000}, {0, True}}], {$p_{(1)}$, 0, 1}, {$p_{(2)}$, 0, $p_{(1)}$}, {$p_{(3)}$, 0, $p_{(2)}$},{$p_{(4)}$, 0, $p_{(3)}$}, {$p_{(5)}$, 0, $p_{(4)}$}, {$p_{(6)}$, 0, $p_{(5)}$}, {$p_{(7)}$, 0, $p_{(6)}$}, {$p_{(8)}$, 0, $p_{(7)}$}, {$p_{(9)}$, 0, $p_{(8)}$}, {$p_{(10)}$, 0, $p_{(9)}$}]

As we see, the above integral is the same as the previous one, exept the multiplier ($p_{(1)}$+ 2*$p_{(2)}$ + 3*$p_{(3)}$ + 4*$p_{(4)}$ + 5*$p_{(5)}$ + 6*$p_{(6)}$ + 7*$p_{(7)}$ + 8*$p_{(8)}$ + 9*$p_{(9)}$ + 10*$p_{(10)}$) of the integrand. Therefore, we need to sum up p1+2*p2+3*p3+4*p4+5*p5+6*p6+7*p7+8*p8+9*p9+10*p10 over all the permutations of indices 1,2,...10. The result is 11!/2*(p1+p2+p3+p4+p5+p6+p7+p8+p9+p10) (consider p1+2*p2+3*p3 as a model). Now we execute the first 7 integrations symbolically

a = Integrate[ 11!/2*(p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9 + p10)*
Piecewise[{{1,  p1*p2 <= 1/400 || ((-p3 + p2) (-p2 + p1))/((1 - p3) (1 - p2)) <= 1/1000000}, {0, True}}],
 {p4, 0, p3}, {p5, 0, p4}, {p6, 0,  p5}, {p7, 0, p6}, {p8, 0, p7}, {p9, 0, p8}, {p10, 0, p9}]

1980*p3^7*(2*p1 + 2*p2 + 9*p3)* Piecewise[{{1, 400*p1*p2 <= 1 || ((p1 - p2)*(p2 - p3))/((-1 + p2)*(-1 + p3)) <= 1/1000000}}, 0]

and the last three integrations numerically

NIntegrate[1980*p3^7*(2*p1 + 2*p2 + 9*p3)*
Piecewise[{{1,  400*p1*p2 <= 1 || ((p1 - p2)*(p2 - p3))/((-1 + p2)*(-1 + p3)) <=  1/1000000}}, 0], {p1, 0, 1}, {p2, 0, p1}, {p3, 0,  p2}]

0.000724406

The same result is produced with the Method->"GlobalAdaptive" option.