Re: Accuracy problem in Mathematica

*To*: mathgroup at smc.vnet.net*Subject*: [mg49118] Re: Accuracy problem in Mathematica*From*: ab_def at prontomail.com (Maxim)*Date*: Fri, 2 Jul 2004 02:01:18 -0400 (EDT)*References*: <cbu21o$5ea$1@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

"aaaa" <aaa at huji.ac.il> wrote in message news:<cbu21o$5ea$1 at smc.vnet.net>... > Hello, > > I'm having a problem with a calculation in Mathematica which I can't > solve. I have an expression which I know to be (from analytical reasons) > always between 0 and 1. It's a function of a and n ( n being natural and > a rational) and it looks like this: > > 1/(1-a^2)^n + > > Sum[((2*n - k - 1)!/((n - 1)!*(n - k)!*2^(2*n - k)))*(1/(1 + a)^k - 1/(1 > - a)^k), {k, 1, n}] > > > > Let's say a=0.5. > > Now, when I try to calculate for small n, it's ok. When calculating for > large n's (around 400 and above) I'm starting to get wrong results (the > number not being between 0 and 1). The problem is that the first term > (the first line before the sum) is VERY VERY close to the negative of > the second term (the sum), and it's getting closer as n grows. When > using large n's, Mathematica says they are the same number or even that > the last term is bigger (which means the whole expression becomes > negative) - which is wrong. It's a matter of accuracy, and I'm not sure > how I can fix it. > > Can anybody help me? > > > > Itamar The straightforward way is to perform all calculations using arbitrary-precision numbers: In[1]:= f[a_, n_]:= 1/(1-a^2)^n + NSum[((2*n - k - 1)!/ ((n - 1)!*(n - k)!*2^(2*n - k)))*(1/(1 + a)^k - 1/(1 - a)^k), {k, 1, n}, Method->Fit, NSumTerms->n, WorkingPrecision->150] f[1/2, 1000] // Timing Out[2]= {7.172 Second, 0.0356247970508144815921869} Actually in this case it wouldn't make much difference if we simply used Sum with high-precision values of a and n. As this would also take significant time, there is a more contrived way which relies on evaluating the sum symbolically: In[3]:= exact = 1/(1-a^2)^n + Sum[ ((2*n - k - 1)!/ ((n - 1)!*(n - k)!*2^(2*n - k)))*(1/(1 + a)^k - 1/(1 - a)^k), {k, 1, n}] Out[3]= (1 - a^2)^(-n) - (2^(1 - 2*n)*Pi*Csc[2*n*Pi]* (HypergeometricPFQRegularized[{1, 1 - n}, {2 - 2*n}, -(2/(-1 + a))] + a*HypergeometricPFQRegularized[{1, 1 - n}, {2 - 2*n}, -(2/(-1 + a))] - HypergeometricPFQRegularized[{1, 1 - n}, {2 - 2*n}, 2/(1 + a)] + a*HypergeometricPFQRegularized[{1, 1 - n}, {2 - 2*n}, 2/(1 + a)]))/ ((-1 + a)*(1 + a)*(-1 + n)!*Gamma[n]) The problem is that Mathematica returns an expression which is undefined for integer values of n, it is correct only as the limit when n->n0 and n0 is integer. What we can do is to evaluate this expression slightly off the integer point: In[4]:= f2 = Function @@ {{a,n}, exact /. n -> n + 1`150*^-150}; Timing[f2[1/2, 1000]] Out[5]= {0.594 Second, 0.0356247970508144815922 + 0.``172*I} Taking the 1`150*^-150 magic number is of course just heuristics, but it seems to work quite well. As a side note, here's why I used the Function@@L method above. Compare the next two outputs: In[6]:= Clear[y, func] y = x; func[x_] = y; DownValues[func] Out[9]= {HoldPattern[func[x_]] :> x} In[10]:= Clear[y, func] Module[{y = x}, func[x_] = y]; DownValues[func] Out[12]= {HoldPattern[func[x$_]] :> x} The resulting functions are different; in my opinion, this is an inconsistency, so I just try to avoid such constructs whenever possible, and Function@@{x,expr} or Set@@{f[x_],expr} bypasses the variable renaming. Maxim Rytin m.r at inbox.ru