Re: Integrating "If"
- To: mathgroup at smc.vnet.net
- Subject: [mg84555] Re: Integrating "If"
- From: "David W.Cantrell" <DWCantrell at sigmaxi.net>
- Date: Thu, 3 Jan 2008 20:29:17 -0500 (EST)
- References: <flc7vh$809$1@smc.vnet.net> <flfubb$jrm$1@smc.vnet.net> <flie4b$fit$1@smc.vnet.net>
hemphill at alumni.caltech.edu wrote:
> "David W.Cantrell" <DWCantrell at sigmaxi.net> writes:
> > hemphill at alumni.caltech.edu wrote:
[snip]
> >> Keep in mind that I am also interested in evaluating integrals like:
> >>
> >> Integrate[If[(x0-x1)^2+(y0-y1)^2<d^2,1,0],{x0,0,1},{x1,0,1},{y0,0,1},
> >> {y1, 0,1}]
> >>
> >> (This integral takes a while to evalulate, and the result isn't
> >> continuous or monotonic. I'm about half way done with evaluating it
> >> by hand--it's an interesting problem.)
> >
> > Yes, it does take a while, and I didn't have the patience to wait. But
> > at least you can see the result fairly quickly using, say,
> >
> > Plot[NIntegrate[UnitStep[d^2 - (x0 - x1)^2 - (y0 - y1)^2], {x0, 0, 1},
> > {x1, 0, 1}, {y0, 0, 1}, {y1, 0, 1}], {d, 0, 2}]
>
> This didn't work for me. The integrator complained that the result
> coverged too slowly.
Perhaps you will be comforted to know that the plot I mentioned above
appears to be the same as a plot obtained from your solution below.
> > BTW, what did you get by hand eventually?
>
> If 0 <= d <= 1, then
>
> (-8*d^3)/3 + d^4/2 + d^2*Pi
>
> If 1 <= d <= Sqrt[2], then
>
> 1/3 - 2*d^2 - d^4/2 + (4*Sqrt[-1 + d^2]*(1 + 2*d^2))/3 + Pi +
> 2*(-1 + d^2)*ArcCot[Sqrt[-1 + d^2]] - 2*(1 + d^2)*ArcTan[Sqrt[-1 +
> d^2]]
>
> I still used Mathematica to do some of the algebra. I had to coerce
> Mathematica to give the answer in this form, since it wanted to return
> answers with "I" and "ArcSinh" in them.
Thanks for showing your result and explaining how you got it.
BTW, in case you're not aware of it, the last three terms of your second
case can be simplified somewhat. Namely, for the values of d under
consideration,
Pi + 2*(-1+d^2)*ArcCot[Sqrt[-1+d^2]] - 2*(1+d^2)*ArcTan[Sqrt[-1+d^2]]
is the same as
d^2 (Pi - 4 ArcTan[Sqrt[d^2 - 1]])
but I know of no direct way to demonstrate that fact using Mathematica.
(Mathematica can show that the derivative of the difference of the two
expressions is 0, etc. But that's hardly a direct method.)
David
> I had started to use a complicated method which evaluated areas of
> circles which were clipped by the unit square. Maxim Rytin (private
> corresponence) pointed out that the pdf of (x0-x1)^2 is (1/Sqrt[u]-1),
> 0<u<=1, and to get the pdf of (x0-x1)^2+(y0-y1)^2, you merely have to
> convolve the first pdf with itself. After that, you integrate to get
> the cdf, and substitute d^2 for u. A little messy, but not too bad.
>
> Scott