[Date Index] [Thread Index] [Author Index]

Re: Bug Report - Two numerical values for a same variable

• To: mathgroup at smc.vnet.net
• Subject: [mg54255] Re: Bug Report - Two numerical values for a same variable
• From: "Carl K. Woll" <carlw at u.washington.edu>
• Date: Mon, 14 Feb 2005 21:50:50 -0500 (EST)
• Organization: University of Washington
• References: <cuhr4q$978$1@smc.vnet.net> <cukc65$lup$1@smc.vnet.net> <200502130521.AAA03676@smc.vnet.net> <opsl4914q7iz9bcq@monster> <002301c51217$ad7438b0$b5803e44@Dell> <cup6bc$e11$1@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

Bob,

Objects with head Rational are atoms, so you don't have access to the 
integers inside a rational object. It's the samething as complex numbers. 
For example,

Complex[6,7] /. 7->1
6 + 7 I

Carl Woll

"DrBob" <drbob at bigfoot.com> wrote in message 
news:cup6bc$e11$1 at smc.vnet.net...
> 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: [mg54255] Re:  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>