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