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