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
    $\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 '19 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$ – Logan Smith Nov 21 '19 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 '19 at 1:08
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}$

| improve this answer | |
$\endgroup$
  • $\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$ – Logan Smith Dec 6 '19 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.