Mathematica 9 is now available
Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2004
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2004

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

Search the Archive

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


  • Prev by Date: Re: Accuracy problem in Mathematica
  • Next by Date: FW: Accuracy problem in Mathematica
  • Previous by thread: Re: Accuracy problem in Mathematica
  • Next by thread: FW: Accuracy problem in Mathematica