FW: Accuracy problem in Mathematica

*To*: mathgroup at smc.vnet.net*Subject*: [mg49125] FW: [mg49061] Accuracy problem in Mathematica*From*: "Wolf, Hartmut" <Hartmut.Wolf at t-systems.com>*Date*: Fri, 2 Jul 2004 02:01:26 -0400 (EDT)*Sender*: owner-wri-mathgroup at wolfram.com

-----Original Message----- From: Wolf, Hartmut To: mathgroup at smc.vnet.net Subject: [mg49125] RE: [mg49061] Accuracy problem in Mathematica >-----Original Message----- >From: aaaa [mailto:aaa at huji.ac.il] To: mathgroup at smc.vnet.net >Sent: Wednesday, June 30, 2004 11:34 AM >To: mathgroup at smc.vnet.net >Subject: [mg49125] [mg49061] Accuracy problem in Mathematica > > >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 > > Itamar, when I evaluate your expression (with n, a undefined) Mathematica simplifies it to In[120]:= expr = 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[120]= (1 - a^2)^(-n) + (1/(-1 + n)!)* ((-((2*Pi*Csc[2*n*Pi]*HypergeometricPFQRegularized[ {1, 1 - n}, {2 - 2*n}, 2/(1 - a)])/ ((-1 + a)*Gamma[n])) - (2*Pi*Csc[2*n*Pi]*HypergeometricPFQRegularized[ {1, 1 - n}, {2 - 2*n}, 2/(1 + a)])/ ((1 + a)*Gamma[n]))/2^(2*n)) Now giving n, a values: In[133]:= n = 20; a = 1/2; expr >From In[133]:= \[Infinity]::"indet": "Indeterminate expression 0 Pi ComplexInfinity encountered." >From In[133]:= \[Infinity]::"indet": "Indeterminate expression 0 Pi ComplexInfinity encountered." Out[133]= Indeterminate This is because the terms: In[134]:= Csc[2 n Pi] Out[134]= ComplexInfinity In[135]:= HypergeometricPFQRegularized[{1, 1 - n}, {2 - 2 n}, 2(1 - a)] Out[135]= 0 So the Simplification of Mathematica is not usable in your case. (This appears to be for all positive Interger n; maybe you're able to proof that) But your expression is finite and rational when n, a get Integer and Rational values before the summation, e.g.: In[136]:= Block[{n = 400, a = 1/2}, 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[136]= 516307148039893425547168143532290205071415791918265630465198976168080842750613\ 809499416727419043626362456623184319361397387028360858469908614633295744633328\ 862859160709473304622880050348452004683108502066519849464065629946092128637050\ 746846813671905860596256092371406885602718086010232583371639746378672299278469\ 876235477885047244577634657496497715524728140641132606408239769827006443281375\ 57151280261324604136115868279422184239/\ 918815807062947172002982616630447315694409668429175232919346394540127635379951\ 214984232652958021474172881030243332247949592807000680057999868079932529077985\ 374292871633644370258227968982931003670645733650408779245993508719468999377001\ 957500228588384000620222104017057282938170614727767933541288546863279277349439\ 745536223805363981702198726988983034665773318277629283750950553689943505174547\ 935790431818538430777692511101882728448 To get an approximate value just do: In[137]:= % // N Out[137]= 0.0561927 But this doesn't work In[140]:= Block[{n = 400, a = .5}, 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[140]=2.076918743413931`*^34 as the definition a = .5 forces the calculation to be done with machine precision, which goes completely astray, on reasons you stated. -- Hartmut Wolf