Re: Re: [Mathematica 5.1] Bug Report - Two numerical values for a same variable
- To: mathgroup at smc.vnet.net
- Subject: [mg54229] Re: [mg54195] Re: [Mathematica 5.1] Bug Report - Two numerical values for a same variable
- From: DrBob <drbob at bigfoot.com>
- Date: Sun, 13 Feb 2005 22:17:19 -0500 (EST)
- References: <cuhr4q$978$1@smc.vnet.net> <cukc65$lup$1@smc.vnet.net> <200502130521.AAA03676@smc.vnet.net> <opsl4914q7iz9bcq@monster> <002301c51217$ad7438b0$b5803e44@Dell>
- Reply-to: drbob at bigfoot.com
- Sender: owner-wri-mathgroup at wolfram.com
I'm not sure why Maxim's method doesn't encounter the bug.
But this seems even more odd, to me:
f[w_] = (Cos[w] - Cos[(23/100)*w]/E^(w^2/2))/w;
a = E^Integrate[f[w], {w, 0, Infinity}];
a /. 20000 -> 20000.`20.
E^((1/2)*(-EulerGamma - Log[2] -
Derivative[1, 0, 0][
Hypergeometric1F1][0, 1/2,
-(529/20000)]))
Cases[a, 20000, Infinity]
{}
FullForm[a]
FullForm[E^(Rational[1, 2]*
(-EulerGamma - Log[2] -
Derivative[1, 0, 0][
Hypergeometric1F1][0,
Rational[1, 2], Rational[
-529, 20000]]))]
If ReplaceAll can find -529/20000 in a, why can't it (or Cases) find 20000?
Bobby
On Sun, 13 Feb 2005 22:02:14 -0000, David W. Cantrell <DWCantrell at sigmaxi.org> wrote:
> ----- Original Message -----
> From: "DrBob" <drbob at bigfoot.com>
To: mathgroup at smc.vnet.net
> To: "David W. Cantrell" <DWCantrell at sigmaxi.org>; <mathgroup at smc.vnet.net>
> Sent: Sunday, February 13, 2005 16:43
> Subject: [mg54229] Re: [mg54195] Re: [Mathematica 5.1] Bug Report - Two numerical values for a same variable
>
>
>> >> But wasn't Ling _already_ using arbitrary-precision numbers in the input?
>>
>> No, he was using exact numbers, and got an answer that was still exact. The substitution
>>
>> -529/20000 -> -529/20000`20
>>
>> made the expression approximate, but more than machine-precision. Numbers like 20000`20 are called "arbitrary precision" for some unfathomable reason.
>
>
> Hmm. Thanks for the info. For some reason (apparently, not a good one), I must have thought that "exact" and "arbitrary-precision" were synonymous.
>
> In any event, this thread has spurred me to upgrade to version 5.1. Having done so, I now see that a bug which I had diagnosed in version 5.0 is still present in version 5.1.
>
> In[1]:=
> Exp[Integrate[(Cos[w] - Exp[-(w^2/2)]*Cos[(23/100)*w])/w, {w, 0, Infinity}]]
> Out[1]=
> E^((1/2)*(-EulerGamma - Log[2] - Derivative[1, 0, 0][Hypergeometric1F1][0, 1/2, -(529/20000)]))
> Version 5.0's exact result was incorrect because it was missing -EulerGamma, BTW.
>
> In[2]:= N[%]
> Out[2]= 0.5439146227050919
>
> which is correct.
>
> In[3]:= N[%%, 20]
> Divide :: infy : Infinite expression 1/0. encountered. ...
> Out[3]=
> 0.52983935469483822112357538077323822033`20.000000000000007
>
> Note that this is incorrect for exactly the same reason which I diagnosed in my original response in this thread: Mathematica treats the derivative of the hypergeometric function as if it were 0. It's easy to see this because Out[6] is the same as Out[3]. So this is a bug present in both versions 5.0 and 5.1.
>
> In[6]:=
> N[E^((1/2)*(-EulerGamma - Log[2])), 20]
> Out[6]=
> 0.52983935469483822112357538077323822033`20.000000000000007
>
> But note that, when asked about that derivative separately, Mathematica knows that it's not 0:
>
> In[8]:=
> N[Derivative[1, 0, 0][Hypergeometric1F1][0, 1/2, -(529/20000)], 20]
> Out[8]=
> -0.05243686946160788253943021003847243652`20.
>
> Finally, I'm rather surprised that Maxim's
>
>> >> f[w_] = (Cos[w] - E^(-w^2/2)*Cos[23/100*w])/w;
>> >> a = E^Integrate[f[w], {w, 0, Infinity}];
>> >> a /. -529/20000 -> -529/20000`20
>
> doesn't suffer from the same bug. Why does his method work correctly, while N[a, 20] encounters the bug?
>
> David Cantrell
>
>
>
>> On Sun, 13 Feb 2005 00:21:26 -0500 (EST), David W. Cantrell <DWCantrell at sigmaxi.org> wrote:
>>
>> > Maxim <ab_def at prontomail.com> wrote:
>> >> On Fri, 11 Feb 2005 08:42:02 +0000 (UTC), Ling Li <ling at caltech.edu>
>> >> wrote:
>> >>
>> >> > Hi,
>> >> >
>> >> > Mathematica v5.1 is much better than v5.0 on the accuracy of
>> >> > integration. I am not talking about numerical accuracy---v5.0 sometimes
>> >> > gives out two answers that are off by a constant for the same input.
>> >> >
>> >> > However, I still get two different numerical values for a same variable
>> >> > in v5.1
>> >> >
>> >> > a=Exp[Integrate[(Cos[w]-Exp[-w^2/2]*Cos[23/100*w])/w, {w,0,Infinity}]]
>> >> > N[a]
>> >> > N[a,20]
>> >> >
>> >> > Thanks for looking into this problem.
>> >> > --Ling
>> >
>> > I apologize for not having checked the result given in my previous post
>> > via numerical integration. Had I done so, I would have realized that the
>> > result was incorrect. Your result of approx 0.5439146 indeed appears to be
>> > correct.
>> >
>> > Problems may be compounded by the fact that I'm still using version 5.0.
>> >
>> >> The expression which is evaluated incorrectly is N[a, 20] (I get
>> >> Divide::infy warnings and 0.52983935469483822112 as the output). This is
>> >> not so unusual when Hypergeometric and MeijerG are involved; for example,
>> >> see http://forums.wolfram.com/mathgroup/archive/2004/Jul/msg00185.html .
>> >> If we start with arbitrary-precision numbers in the input, we obtain the
>> >> correct result without any problems:
>> >
>> > But wasn't Ling _already_ using arbitrary-precision numbers in the input?
>> >
>> > David
>> >
>> >> In[1]:=
>> >> f[w_] = (Cos[w] - E^(-w^2/2)*Cos[23/100*w])/w;
>> >> a = E^Integrate[f[w], {w, 0, Infinity}];
>> >> a /. -529/20000 -> -529/20000`20
>> >>
>> >> Out[3]=
>> >> 0.54391462270509428
>> >>
>> >> Here are two ways to verify the result. First, you can use the lim
>> >> function from
>> >> http://forums.wolfram.com/mathgroup/archive/2003/Aug/msg00233.html to
>> >> simplify the derivative of Hypergeometric1F1:
>> >>
>> >> In[7]:=
>> >> a /. Derivative[1, 0, 0][Hypergeometric1F1][a_, b_, z_] :>
>> >> lim[(Hypergeometric1F1[a + eps, b, z] -
>> >> Hypergeometric1F1[a, b, z])/eps,
>> >> eps -> 0]
>> >> N[%, 20]
>> >>
>> >> Out[7]=
>> >> E^((-EulerGamma + (529*HypergeometricPFQ[{1, 1}, {3/2, 2},
>> >> -529/20000])/10000 - Log[2])/2)
>> >>
>> >> Out[8]=
>> >> 0.54391462270509428214
>> >>
>> >> Another way is to evaluate the integral numerically (well, almost). To do
>> >> that, we will split the integration range into [0, 1] and [1, Infinity)
>> >> and then integrate the term with E^(-w^2/2) using Method ->
>> >> DoubleExponential and integrate Cos[w]/w symbolically:
>> >>
>> >> In[9]:=
>> >> E^(NIntegrate[f[w], {w, 0, 1}, WorkingPrecision -> 30] +
>> >> (If[FreeQ[#, E],
>> >> Integrate[#, {w, 1, Infinity}],
>> >> NIntegrate[#, {w, 1, Infinity},
>> >> Method -> DoubleExponential, WorkingPrecision -> 30]
>> >> ]& /@ Expand[f[w]]))
>> >>
>> >> Out[9]=
>> >> 0.5439146227050942821365022967
>> >>
>> >> We could also do it in a purely numerical way, integrating the
>> >> oscillatory terms on [1, Infinity) separately by using Method ->
>> >> Oscillatory and increasing the value of MaxRecursion, but I wasn't able
>> >> to get the result with more than 15 digits of precision then.
>> >>
>> >> Mathematica has trouble with high-precision evaluation of MeijerG as
>> >> well: N[MeijerG[{{0, 1}, {}}, {{1/2}, {}}, 1/100], 20] seems unable to
>> >> find even one significant digit with any setting of $MaxExtraPrecision,
>> >> and gives an incorrect result if we increase $MinPrecision (MeijerG[{{0,
>> >> 1}, {}}, {{1/2}, {}}, z] is O[Sqrt[z]] when z->0). Perhaps there are no
>> >> general enough methods to evaluate MeijerG in such cases, but the claim
>> >> that Mathematica "can evaluate special functions with any parameters to
>> >> any precision" seems a little too strong anyway.
>
--
DrBob at bigfoot.com
www.eclecticdreams.net
- References:
- Re: [Mathematica 5.1] Bug Report - Two numerical values for a same variable
- From: "David W. Cantrell" <DWCantrell@sigmaxi.org>
- Re: [Mathematica 5.1] Bug Report - Two numerical values for a same variable