Re: Multiple integration of UnitStep fails
- To: mathgroup at smc.vnet.net
- Subject: [mg63321] Re: Multiple integration of UnitStep fails
- From: Maxim <m.r at inbox.ru>
- Date: Fri, 23 Dec 2005 05:08:42 -0500 (EST)
- References: <do0n3j$4ku$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
On Sat, 17 Dec 2005 09:51:47 +0000 (UTC), Dr. Wolfgang Hintze
<weh at snafu.de> wrote:
> Hello group,
>
> trying to solve the nice problem of determining the probability pn that
> a polygon formed by n (>=4) random points on the unit circle is void of
> an acute angle I came up with the following multiple integral (written
> down here for n=5)
>
> In[15]:=
> p5 = (4!*Integrate[Integrate[Integrate[Integrate[UnitStep[\[Phi]4 - Pi,
> \[Phi]5 - \[Phi]2 - Pi, Pi - \[Phi]3, Pi - \[Phi]4 + \[Phi]2,
> Pi - \[Phi]5 + \[Phi]3], {\[Phi]5, \[Phi]4, 2*Pi}], {\[Phi]4,
> \[Phi]3, 2*Pi}], {\[Phi]3, \[Phi]2, 2*Pi}], {\[Phi]2, 0, 2*Pi}])/(2*Pi)^4
>
> Mathematica version 4 was not able to solve this but returned it
> unevaluated after some minutes; version 5 complained several things
> like: argument is not a power series, unable to check convergence, but
> didn't come up with any result in ten minutes (I wouldn't wait longer).
>
> I could successfully check the normalization at least:
>
> In[14]:=
> p5 = (4!*Integrate[Integrate[Integrate[Integrate[UnitStep[1], {\[Phi]5,
> \[Phi]4, 2*Pi}], {\[Phi]4, \[Phi]3, 2*Pi}], {\[Phi]3, \[Phi]2, 2*Pi}],
> {\[Phi]2, 0, 2*Pi}])/(2*Pi)^4
>
> Out[14]= 1
>
> How would you proceed to solve In[15]? What about the general case (n=6
> see below)?
>
> Any hints are greatly appreciated.
>
> Regards,
> Wolfgang
>
> PS:
>
> Here's the probability for the case n=6
>
> p6 =
> (5!*Integrate[Integrate[Integrate[Integrate[Integrate[UnitStep[\[Phi]5 -
> Pi, \[Phi]6 - \[Phi]2 - Pi, Pi - \[Phi]3, Pi - \[Phi]4 + \[Phi]2,
> Pi - \[Phi]5 + \[Phi]3, Pi - \[Phi]6 + \[Phi]4], {\[Phi]6,
> \[Phi]5, 2*Pi}], {\[Phi]5, \[Phi]4, 2*Pi}], {\[Phi]4, \[Phi]3, 2*Pi}],
> {\[Phi]3, \[Phi]2, 2*Pi}], {\[Phi]2, 0, 2*Pi}])/(2*Pi)^5
>
> $Aborted
>
The integral in p6 can be evaluated if we first perform the cylindrical
decomposition of the region where UnitStep is non-zero:
In[1]:= p6 = Hold[
(5!*Integrate[Integrate[Integrate[Integrate[Integrate[
UnitStep[\[Phi]5 - Pi, \[Phi]6 - \[Phi]2 - Pi,
Pi - \[Phi]3, Pi - \[Phi]4 + \[Phi]2,
Pi - \[Phi]5 + \[Phi]3, Pi - \[Phi]6 + \[Phi]4],
{\[Phi]6, \[Phi]5, 2*Pi}],
{\[Phi]5, \[Phi]4, 2*Pi}],
{\[Phi]4, \[Phi]3, 2*Pi}],
{\[Phi]3, \[Phi]2, 2*Pi}],
{\[Phi]2, 0, 2*Pi}])/(2*Pi)^5];
p6 //. HoldPattern[Integrate[Integrate[e_, Siter2__], Siter__]] :>
Integrate[e, Siter, Siter2] /.
HoldPattern[Integrate[UnitStep[s__], Siter__]] :>
Integrate[Boole[Reduce[
Join[(# >= 0&) /@ {s},
(#[[2]] <= #[[1]] <= #[[3]]&) /@ {Siter}],
{Siter}[[All, 1]], Reals]],
Sequence @@ ({#[[1]], -Infinity, Infinity}&) /@ {Siter}] //
ReleaseHold
Out[2]= 1/4
With PiecewiseIntegrate
( http://library.wolfram.com/infocenter/MathSource/5117/ ) we only need to
carry out the first step:
In[3]:= p6 //.
HoldPattern[Integrate[Integrate[e_, Siter2__], Siter__]] :>
Integrate[e, Siter, Siter2] /.
Integrate -> PiecewiseIntegrate //
ReleaseHold
Out[3]= 1/4
Another thing that often works for piecewise integrands is to restrict the
domain of the parameter values:
In[4]:= (int = Integrate[Min[x, y], {x, 0, a}, {y, 0, b}];) // Timing
Out[4]= {114.547*Second, Null}
In[5]:= Cases[int, If[a_, __] :> a, {0, -1}]
Out[5]= {b > 0 && (b <= Re[a] || Re[a] <= 0 || Im[a] != 0)}
In[6]:= int /. {{a -> 2, b -> 1}, {a -> 2, b -> 3}}
Out[6]= {Indeterminate, 6}
We can see that the condition in the answer allows for the possibility of
complex values (for a but not b), which I don't think is in any way
useful, and even at that the conditional answer is incorrect for {a, b} =
{2, 1}. Besides, the condition doesn't cover all the possibilities, so for
{a, b} = {2, 3} the integration has to be performed again, and the result
is also incorrect. With the assumption Element[{a, b}, Reals] the
integration takes significantly less time and we get an answer that is
correct for all values of the parameters.
Maxim Rytin
m.r at inbox.ru