Re: Incorrect integral
- To: mathgroup at smc.vnet.net
- Subject: [mg43642] Re: [mg43623] Incorrect integral
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Sat, 27 Sep 2003 04:58:04 -0400 (EDT)
- References: <40700B4C02ABD5119F0000902787664407D78C09@hplex1.hpl.hp.com>
- Sender: owner-wri-mathgroup at wolfram.com
"Rasmusson, Lars" wrote:
>
> Thank you for help. It seems to work.
>
> However, Mathematica does not seem to be working correctly when the integration boundaries are given symbolicly. Is there a way to handle that as well?
>
> To see that it is incorrect, compare the result of
>
> In[1]:=
> InputForm[f[p_] = Integrate[UnitStep[k-a*p]*(k-a*p),
> {k,K0,K1}, {a,A0,A1}, Assumptions->0<p<1]] ;
>
> In[2]:=
> f[0.1]/.{K0|A0->0,K1|A1->1}
>
> Out[2]=
> 0.453333
>
> with the correct value 0.451667.
>
> Again, I'm very grateful for your good response to my question.
>
> Best regards,
> Lars Rasmusson
>
> -----Original Message-----
> From: Daniel Lichtblau [mailto:danl at wolfram.com]
To: mathgroup at smc.vnet.net
> Sent: Friday, September 26, 2003 7:47 AM
> To: Lars Rasmusson
> Cc: mathgroup at smc.vnet.net
> Subject: [mg43642] Re: [mg43623] Incorrect integral
>
> Lars Rasmusson wrote:
> >
> > Hi, it seems like Integrate is not handling the Max[ ] function
> > properly, or am I mistaken? Compare the two outputs
> > In[1]:=
> >
> > fa[p_] := NIntegrate[Max[0, k - a p], {k, 0, 1}, {a, 0, 1}] fb[p_] :=
> > Integrate[Max[0, k - a p], {k, 0, 1}, {a, 0, 1}] fa[0.1]
> > fb[0.1]
> >
> > Out[3]=
> > 0.451667
> > Out[4]=
> > 1.66667
> >
> > The outputs differ. NIntegrate returns the correct value.
> >
> > In[5]:= $Version
> > Out[5]= 5.0 for Microsoft Windows (June 11, 2003)
> >
> > Is there any way to get a correct behavior from Mathematica? I would
> > be happy to receive the symbolic solution (to a more complicated
> > integral of this kind).
> >
> > Thanks,
> >
> > Lars
>
> One method that Integrate can handle is to reformulate using UnitStep instead of Max.
>
> In[9]:= InputForm[f[p_] = Integrate[UnitStep[k-a*p]*(k-a*p),
> {k,0,1}, {a,0,1}, Assumptions->0<p<1]]
> Out[9]//InputForm= (3 - 3*p + p^2)/6
>
> In[11]:= f[.1]
> Out[11]= 0.451667
>
> Daniel Lichtblau
> Wolfram Research
This one is not so obvious. It has to do with the fact that you get
different limiting behavior depending on how you arrive at the specified
values. In particular we will at first refrain from performing the
substitution for A0.
In[30]:= InputForm[f[p] /. {A1|K1->1, K0->0, p->1/10}]
Out[30]//InputForm=
(5*((27*(-1 + A0/10)*(-1 + A0)*UnitStep[1 - A0/10]*
(UnitStep[1 - A0] + UnitStep[-1 + A0]))/100 +
UnitStep[1 - A0/10]*(((1 - 3*A0 + 3*A0^2)*UnitStep[-A0/10])/1000 -
((-1 + A0)^3*UnitStep[1 - A0]*UnitStep[A0/10])/1000)))/3
Now we do the substitution.
In[4]:= ee /. A0->0
34
Out[4]= --
75
But if you look at the prior expression you see a UnitStep[A0] term
which indicates that we have discontinuous behavior at A0==0. We'll take
a limit as we approach from the positive side (default behavior for
Limit).
In[5]:= Limit[ee, A0->0]
271
Out[5]= ---
600
So you have a path-dependent discontinuity at the intended values. Blind
substitution cannot reliably work in this situation. Offhand I do not
know how to work around this. The discontinuous behavior pretty much
implies that you have to do things carefully at that point.
Daniel Lichtblau
Wolfram Research