Re: Re: [Mathematica 5.1] Bug Report - Two numerical values for a same variable
- To: mathgroup at smc.vnet.net
- Subject: [mg54230] Re: [mg54195] Re: [Mathematica 5.1] Bug Report - Two numerical values for a same variable
- From: "David W. Cantrell" <DWCantrell at sigmaxi.org>
- Date: Sun, 13 Feb 2005 22:17:20 -0500 (EST)
- Sender: owner-wri-mathgroup at wolfram.com
----- Original Message ----- From: "DrBob" <drbob at bigfoot.com> To: mathgroup at smc.vnet.net Subject: [mg54230] 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 reaso= n. 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.