Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

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