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: [mg49117] Re: Accuracy problem in Mathematica
  • From: Paul Abbott <paul at physics.uwa.edu.au>
  • Date: Fri, 2 Jul 2004 02:01:17 -0400 (EDT)
  • Organization: The University of Western Australia
  • References: <cbu21o$5ea$1@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

In article <cbu21o$5ea$1 at smc.vnet.net>, "aaaa" <aaa at huji.ac.il> wrote:

> 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.

One solution is to use high-precision arithmetic. For example, with

  f[n_][a_] := Sum[((-k + 2 n - 1)! (1/(a + 1)^k - 1/(1 - a)^k))/
    ((n - 1)! (n - k)! 2^(2 n - k)), {k, 1, n}] + 1/(1 - a^2)^n

then

  f[400][0.5`70.]

is between 0 and 1 and gives you a precision of just over 16.

Note that you can change the upper limit of the summation to Infinity 
(because of the (n - k)! term in the denominator) and obtain an 
alternative form for the sum:

  g[n_][a_] = Simplify[Sum[((-k + 2 n - 1)! (1/(a + 1)^k - 1/(1 - a)^k))/
    ((n - 1)! (n - k)! 2^(2 n - k)), {k, 1, Infinity}] + 1/(1 - a^2)^n,
      n \[Element] Integers]

However, high-precision arithmetic is still required:

  g[400][0.5`70.]

The optimal solution for calculations with large n would be to determine 
asymptotic expansions of the sum. For example, here is a very crude 
analysis:

  summand[k_] = Simplify[Normal[Series[(2n-k-1)!/
   ((n-1)!(n-k)! 2^(2n - k)), {n, Infinity, 2}]]]

(ignoring the warning messages),

  asympt[n_][a_] = Simplify[Sum[(1/(a + 1)^k - 1/(1 - a)^k) summand[k],
   {k, 1, Infinity}]]

Then compare the sum to the asymptotic expression, computed at low 
precision: 

  g[400][0.5`70.] - asympt[400][0.5]

Cheers,
Paul

-- 
Paul Abbott                                   Phone: +61 8 9380 2734
School of Physics, M013                         Fax: +61 8 9380 1014
The University of Western Australia      (CRICOS Provider No 00126G)         
35 Stirling Highway
Crawley WA 6009                      mailto:paul at physics.uwa.edu.au 
AUSTRALIA                            http://physics.uwa.edu.au/~paul


  • Prev by Date: Re: converting table output and plotting
  • Next by Date: Re: Accuracy problem in Mathematica
  • Previous by thread: Re: Accuracy problem in Mathematica
  • Next by thread: Re: Accuracy problem in Mathematica