MathGroup Archive 2005

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

Search the Archive

Re: [Mathematica 5.1] Bug Report - Two numerical values for a same variable

On Fri, 11 Feb 2005 08:42:02 +0000 (UTC), Ling Li <ling at> 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

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 .  
If we start with arbitrary-precision numbers in the input, we obtain the  
correct result without any problems:

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


Here are two ways to verify the result. First, you can use the lim  
function from to  
simplify the derivative of Hypergeometric1F1:

a /. Derivative[1, 0, 0][Hypergeometric1F1][a_, b_, z_] :>
   lim[(Hypergeometric1F1[a + eps, b, z] -
         Hypergeometric1F1[a, b, z])/eps,
     eps -> 0]
N[%, 20]

E^((-EulerGamma + (529*HypergeometricPFQ[{1, 1}, {3/2, 2},  
-529/20000])/10000 - Log[2])/2)


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:

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]]))


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.

Maxim Rytin
m.r at

  • Prev by Date: Re: Re: how to have a blind factorization of a polinomial?
  • Next by Date: Re: Fourier Transfer and a game?!?!
  • Previous by thread: Re: [Mathematica 5.1] Bug Report - Two numerical values for a same variable
  • Next by thread: Re: [Mathematica 5.1] Bug Report - Two numerical values for a same variable