MathGroup Archive 2005

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

Search the Archive

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.


  • Prev by Date: Re: Re: [Mathematica 5.1] Bug Report - Two numerical values for a same variable
  • Next by Date: Re: rules and lists
  • Previous by thread: Re: Re: Bug Report - Two numerical values for a same variable
  • Next by thread: VertexStyle determined by List elements