6
$\begingroup$

Similar to my previous question, running

Clear["Global`*"];
ineq = (1+y+y^2) r[4]<r[2]+y (y r[1]+r[3]) && y (r[1]-r[2]-(1+y) (r[3]-r[4]))+r[4]<r[2];
reg1=ImplicitRegion[ineq,{{r[1],0,1},{r[2],0,1},{r[3],0,1},{r[4],0,1}}];
(int[y_]=Assuming[1>y>0,Integrate[1,{r[1],r[2],r[3],r[4]}\[Element]reg1]//Simplify])//AbsoluteTiming

Produces $$ \left\{165.709,\frac{8 y^4+26 y^3+39 y^2+32 y+12}{24 (y+1)^2 \left(y^2+y+1\right)}\right\} $$

Is there a way to speed this up using Implicit Region? Or is there a faster method to find integration bounds from an inequality? An explanation for why this isn't possible would also be accepted.

$\endgroup$
3
  • 4
    $\begingroup$ Integration bounds may be obtained from Reduce[1 > y > 0 && First@reg1, {r[1], r[2], r[3], r[4]}] /. _Equal -> False, where the _Equal -> False removes components of measure zero. It will also show some of the complexity. The order of the variable can matter, since you're setting up an iterated integral. $\endgroup$
    – Michael E2
    Nov 20, 2019 at 23:34
  • 1
    $\begingroup$ Doing that shows a hundred lines of complicated integration bounds. Is that the best we can do for this inequality? $\endgroup$ Nov 21, 2019 at 0:43
  • $\begingroup$ There's an unconventional kluge that gets the result a couple of orders of magnitude faster, happy to post if you want, but it is a kluge that is not very generic... $\endgroup$
    – ciao
    Nov 23, 2019 at 1:08

1 Answer 1

7
+50
$\begingroup$

Using @Michael's suggestion brings the timing down about an order of magnitude. There are 120 components after using LogicalExpand:

components = List @@ LogicalExpand[
    Reduce[1 > y > 0 && First @ reg1, {r[1], r[2], r[3], r[4]}] /. _Equal -> False
]; //AbsoluteTiming
Length[components]

{6.53425, Null}

120

Computing integral over each region:

r1 = Assuming[
    0<y<1,
    Simplify @ Total[
        Map[
            Integrate[1, z ∈ ImplicitRegion[#, {r[1], r[2], r[3], r[4]}]]&,
            components
        ]
    ]
]; //AbsoluteTiming

{7.24909, Null}

The result is basically the same as yours, except for a discrete set of values fory:

r1 /. s_Root :> N[s] //TeXForm

$\begin{cases} \frac{8 y^4+26 y^3+39 y^2+32 y+12}{24 (y+1)^2 \left(y^2+y+1\right)} & \frac{1}{2} \left(\sqrt{5}-1\right)<y<0.657298\lor 0.543689<y<\frac{1}{2} \left(\sqrt{5}-1\right)\lor 0.657298<y<0.682328\lor 0.682328<y<0.754878\lor y<0.543689\lor y>0.754878 \\ 0 & \text{True} \end{cases}$

$\endgroup$
1
  • $\begingroup$ This is great for this specific problem, but it fails to generalize. For ineq = r[1] > r[2] , it gives infinity as opposed to 1/2 (half of a square with r[] as the axis). I'm currently trying to understand why it works here, but not in the simpler case. $\endgroup$ Dec 6, 2019 at 17:38

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.