MathGroup Archive 2005

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

Search the Archive

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


  • Prev by Date: Re: Using Select with arrays? (Relative newbie)
  • Next by Date: Re: Bug Report - Two numerical values for a same variable
  • Previous by thread: Re: Axes in ShowGraph
  • Next by thread: Re: Bug Report - Two numerical values for a same variable