Re: Wrong Integral result for a Piecewise function
- To: mathgroup at smc.vnet.net
- Subject: [mg58561] Re: Wrong Integral result for a Piecewise function
- From: Maxim <ab_def at prontomail.com>
- Date: Thu, 7 Jul 2005 05:35:51 -0400 (EDT)
- References: <dag1bg$5j4$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
On Wed, 6 Jul 2005 07:30:56 +0000 (UTC), Dean Nairn <dnairn at udel.edu>
wrote:
> This integral gives gives the wrong result for the interval [2,3]
>
> h[x_] := Integrate[Boole[x - 1 < 2 y + 2 z < x], {y, 0, 1}, {z, 0, 1}]
>
> Plot[h[x],{x,0,5}]
> and
> Plot[Evaluate[h[x]],{x,0,5}]
>
> give different plots, The curve should be smooth and bell shaped from
> 0 to 5, so the first looks correct. The second has a jump discontinuity
> at 2 and 3. Also
>
> h[5/2]
> and
> h[x]/.x-> 5/2
>
> give different answers
>
> This is using some new features in Mathematica 5.1. Same result on a
> Mac (10.4) and SunOS (5.9).
>
> Breaking into a difference of two integrals gives the correct answer:
> Integrate[Boole[ 2 y + 2 z < x], {y, 0, 1}, {z, 0, 1}] -
> Integrate[Boole[ 2 y + 2 z < x-1], {y, 0, 1}, {z, 0, 1}]
>
> Finally
> Integrate[h[x], {x, 0, 5}] and NIntegrate[h[x], {x, 0, 5}]
> both give the wrong answer, it should be 1. The triple integral is
> correct
>
> Integrate[Boole[x - 1 < 2 y + 2z < x], {y, 0, 1}, {z, 0, 1}, {x, 0, 5}]
>
> Any suggestions on integrating over regions with linear constraints?
> Versions 5.1 has powerful new piecewise functions
>
Have a look at this thread from a few days ago:
http://forums.wolfram.com/mathgroup/archive/2005/Jul/msg00042.html . The
same two recipes apply here: either use PiecewiseIntegrate
In[1]:=
PiecewiseIntegrate[Boole[x - 1 < 2*y + 2*z < x], {y, 0, 1}, {z, 0, 1}]
Out[1]=
If[x == 3, 3/8, 0] + If[0 < x < 1, x^2/8, 0] + If[1 <= x <= 2, (1/8)*(-1
+ 2*x), 0] + If[2 < x < 3, (1/8)*(-9 + 10*x - 2*x^2), 0]
+ If[Inequality[3, Less, x, LessEqual, 4], (1/8)*(9 - 2*x), 0] + If[4 < x
< 5, (1/8)*(-5 + x)^2, 0]
or integrate separately for ranges of the parameter values where the
answer can be given as a single expression without conditions:
In[2]:=
Piecewise[
{Assuming[#,
Integrate[Boole[x - 1 < 2*y + 2*z < x], {y, 0, 1}, {z, 0, 1}]],
#}& /@
(Flatten[#, 1]&)@
Table[{i < x < i + 1, x == i + 1}, {i, 0, 4}]
] // PiecewiseExpand
Out[2]=
Piecewise[{{1/8, x == 1 || x == 4}, {3/8, x == 2 || x == 3}, {(1/8)*(9 -
2*x), 3 < x < 4}, {(1/8)*(-5 + x)^2, 4 < x < 5}, {x^2/8, 0 < x < 1},
{(1/8)*(-1 + 2*x), 1 < x < 2}, {(1/8)*(-9 + 10*x - 2*x^2), 2 < x < 3}}]
This integral also can be derived from the general convolution formula. If
xi_1, ..., xi_n are independent and uniformly distributed on [0, 1], then
the CDF of their sum is given as
F[n_, x_] = (1/n!)*Sum[(-1)^k*Binomial[n, k]*(x - k)^n*UnitStep[x - k],
{k, 0, n}]
Then the integral in question equals F[2, x/2] - F[2, (x - 1)/2].
Also note that if h[x] is called inside Plot, it still uses Integrate but
with an approximate integrand (after substituting a numerical value for
x), which is not the same as using NIntegrate. Besides, NIntegrate in
versions 5.0 and 5.1 pre-evaluates the integrand, so the numerical
integration proceeds with the wrong symbolic value of h[x].
Maxim Rytin
m.r at inbox.ru