Re: Incorrect integral
- To: mathgroup at smc.vnet.net
- Subject: [mg43656] Re: Incorrect integral
- From: "David W. Cantrell" <DWCantrell at sigmaxi.org>
- Date: Sun, 28 Sep 2003 06:00:36 -0400 (EDT)
- References: <40700B4C02ABD5119F0000902787664407D78C09@hplex1.hpl.hp.com> <bl3k6d$eo2$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Daniel Lichtblau <danl at wolfram.com> wrote:
> "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: [mg43656] Re: 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.
Perhaps using UnitStep is the best way currently available to have
Mathematica handle the problem. However the integrand is continuous, and
so it's unfortunate and "artificial" to introduce functions which are
discontinuous.
There is, I think, a result which will allow "blind substitution" to work
reliably in this situation:
(A1 - A0)(K1 - K0)(K0 + K1 - p(A0 + A1))/4 +
(Abs[K1 - p A0]^3 - Abs[K1 - p A1]^3 + Abs[K0 - p A1]^3 - Abs[K0 - p A0]^3)/(12p)
Note, however, that I did not obtain this result using Mathematica.
David Cantrell