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

• To: mathgroup at smc.vnet.net
• Subject: [mg54190] Re: [Mathematica 5.1] Bug Report - Two numerical values for a same variable
• From: Maxim <ab_def at prontomail.com>
• Date: Sat, 12 Feb 2005 01:59:18 -0500 (EST)
• References: <cuhr4q\$978\$1@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

```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
>

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:

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.

Maxim Rytin
m.r at inbox.ru

```

• 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