Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2005
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2005

[Date Index] [Thread Index] [Author Index]

Search the Archive

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


  • Prev by Date: Re: Returning an empty sequence
  • Next by Date: Re: Gaussian sums (Was: Speeding up simple Mathematica expressions?)
  • Previous by thread: Re: Multiple integration of UnitStep fails
  • Next by thread: Re: EUREKA Re: Types in Mathematica, a practical example