Re: Bug Report - Two numerical values for a same variable
- To: mathgroup at smc.vnet.net
- Subject: [mg54261] Re: Bug Report - Two numerical values for a same variable
- From: DrBob <drbob at bigfoot.com>
- Date: Mon, 14 Feb 2005 21:50:59 -0500 (EST)
- References: <00ed01c512b0$2f242850$6400a8c0@Main>
- Reply-to: drbob at bigfoot.com
- Sender: owner-wri-mathgroup at wolfram.com
That explains it, but only in the sense that "things fall down" is a theory of gravity. Why should Rationals be atomic, for goodness sake? And how did I use Mathematica all this time without hearing about it? Sigh... Bobby On Mon, 14 Feb 2005 11:13:53 -0500, Carl K. Woll <carlw at u.washington.edu> wrote: > 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: [mg54261] 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 >> > > > > > > > -- DrBob at bigfoot.com www.eclecticdreams.net