Re: Accuracy problem in Mathematica

*To*: mathgroup at smc.vnet.net*Subject*: [mg49108] Re: [mg49061] Accuracy problem in Mathematica*From*: Garry Helzer <gah at math.umd.edu>*Date*: Thu, 1 Jul 2004 05:26:06 -0400 (EDT)*References*: <200406300934.FAA05317@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

On Jun 30, 2004, at 5:34 AM, aaaa wrote: > 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 > > Well, you can increase accuracy. In fact you can calculate it exactly and then convert to floating point. g[a_, n_] := 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}] In[25]:= g[5/10,400] Out[25]= 516307148039893425547168143532290205071415791918265630465198976168080842 750613\ 809499416727419043626362456623184319361397387028360858469908614633295744 633328\ 862859160709473304622880050348452004683108502066519849464065629946092128 637050\ 746846813671905860596256092371406885602718086010232583371639746378672299 278469\ 876235477885047244577634657496497715524728140641132606408239769827006443 281375\ 57151280261324604136115868279422184239/\ 918815807062947172002982616630447315694409668429175232919346394540127635 379951\ 214984232652958021474172881030243332247949592807000680057999868079932529 077985\ 374292871633644370258227968982931003670645733650408779245993508719468999 377001\ 957500228588384000620222104017057282938170614727767933541288546863279277 349439\ 745536223805363981702198726988983034665773318277629283750950553689943505 174547\ 935790431818538430777692511101882728448 In[26]:= %//N Out[26]= 0.0561927 Interestingly, Mathematica 5.0 gives a closed form for this function In[28]:= f[a_,n_]=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[28]= \!\(\((1 - a\^2)\)\^\(-n\) - \((2\^\(1 - 2\ n\)\ ¹\ Csc[2\ n\ ¹]\ \((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])\)\) but it has problems: In[27]:= f[5/10,400] From In[27]:= \!\(\* RowBox[{\(°::"indet"\), \(\(:\)\(\ \)\), "\<\"Indeterminate expression \\!\\(0\\\\ ¹\\\\ ComplexInfinity\\) encountered. \\!\\(\\*ButtonBox[\\\"MoreÉ\\\", \ ButtonStyle->\\\"RefGuideLinkText\\\", ButtonFrame->None, \ ButtonData:>\\\"General::indet\\\"]\\)\"\>"}]\) Out[27]= Indeterminate Garry Helzer gah at math.umd.edu