Re: Numerical accuracy of Hypergeometric2F1
- To: mathgroup at smc.vnet.net
- Subject: [mg55862] Re: Numerical accuracy of Hypergeometric2F1
- From: "Christos Argyropoulos M.D." <chrisarg at fuse.net>
- Date: Fri, 8 Apr 2005 01:36:37 -0400 (EDT)
- References: <d3043c$nei$1@smc.vnet.net> <paul-63DCFE.23493207042005@news.uwa.edu.au>
- Sender: owner-wri-mathgroup at wolfram.com
Thanx for the input.
The reason i repeated the "experiment" was to understand what exactly is
going on, when Mathematica evaluates the expression for the hypergeometric
so I can betchmark a Non-Mathematica application using this function .
Thanks also for pointing out that when I asked for the answer to
Hypergeometric2F1[N[90.+1/2,1000],1,3/2,N[90/100,1000], I was not really
getting "an arbitrary precision" approximation to the hypergeometric
function.
Using
Hypergeometric2F1[SetPrecision[90.+1/2,1000],1,3/2,SetPrecision[90/100,1000]]
gives the "correct answer".
Christos Argyropoulos
----- Original Message -----
From: "Paul Abbott" <paul at physics.uwa.edu.au>
To: mathgroup at smc.vnet.net
Subject: [mg55862] Re: Numerical accuracy of Hypergeometric2F1
> In article <d3043c$nei$1 at smc.vnet.net>,
> "Christos Argyropoulos M.D." <chrisarg at fuse.net> wrote:
>
>> Thanx everyone for pointing out the subtleties of making Mathematica
>> understand precision of input (I always thought that SetPrecision threads
>> over its first input).
>
> Why?
>
>> However I still think that the asymptotic expansion
>> (or whatever numerical method is used internally) used to obtain the
>> numerical value of Hypergeometric2F1 for large (floating point, arbitrary
>> precision) values of one of the parameters needs to be re-examined by
>> WRI.
>> I am attaching an example where the parameters to hypergeometric2F1 were
>> arbitrary precision numbers (using SetPrecision before passing them to
>> the
>> function).
>
> In that example you set
>
> k=90.;
>
> This is a machine-precision number. Using it in artihmetic with other
> expressions will coerce them to machine-precision also. For example, try
>
> k Pi
>
> or
>
> 2 k
>
> or
>
> 2.7`30 k
>
> and you will get back a machine-precision number -- so your example does
> not quite work the way you think that it does.
>
>> When rational numbers or integers are used to represent input,
>> Mathematica
>> gets it right (I believe that in those cases it symbolically expands the
>> hypergeometric in a polynomial and then evaluates the polynomial).
>> For the example attached (i.e. Hypergeometric2F1[80+1/2,1,3/2,90/10])
>
> You mean Hypergeometric2F1[90+1/2,1,3/2,90/10]. If you enter
>
> Hypergeometric2F1[90+1/2,1,3/2,x]
>
> you will see that Mathematica indeed expands the hypergeometric into a
> polynomial.
>
>> and floating point representations for [90+1/2, 90/10] one has to
>> increase
>> precision to 25 digits to get an accurate result (and relying on machine
>> arithmetic is off the mark). By comparison, Forrey's implementation of
>> the
>> Hypergeometric2F1 ( R. C. Forrey, Computing the Hypergeometric Function,
>> Journal of Computational Physics 137, 79-100 (1997). ) which uses machine
>> floating point arithmetic gives the following answer (also obtainable by
>> Mathematica for integer/rational input) :
>>
>> 0.9860651611217461E+89
>
> Actually, Mathematica gives
>
> N[Hypergeometric2F1[90 + 1/2, 1, 3/2, 9/10]]//InputForm
> 9.860651611217265*^88
>
> and computation to higher precision confirms that all digits except the
> last are correct -- so Forrey's result is also correct to only 13
> digits.
>
>> Using a smaller first argument (i.e. 10) results in no "significant"
>> difference in the answers provided by Mathematica irrespective of the
>> encoding of the parameters (rational, machine arithmetic floating point,
>> arbitrary precision numbers). My interpretation of these results is that
>> the
>> algorithm that provides the numerical answer needs some
>> tweaking/optimization
>
> The linear transformation used for 1/2 < x <= 1 is
>
> Hypergeometric2F1[k + 1/2, 1, 3/2, x] /.
> Hypergeometric2F1[a_, b_, c_, z_] :>
> (Gamma[c] Gamma[c - a - b])/(Gamma[c - a] Gamma[c -b]) *
> Hypergeometric2F1[a, b, a + b - c + 1, 1 - z] +
> (Gamma[c] Gamma[a + b - c])/(Gamma[a] Gamma[b]) *
> Hypergeometric2F1[c - a, c - b, c - a - b + 1, 1 - z]/
> (1 - z)^(a + b - c) // FullSimplify
>
> and, for large k, using only the term
>
> largek = (Sqrt[Pi] Gamma[k])/((1 - x)^k 2 Sqrt[x] Gamma[k + 1/2])
>
> with machine precision values also gives 13 correct digits
>
> largek /. {x -> 0.9, k -> 90.} // InputForm
>
> 9.860651611217077*^88
>
> (the second term is insignificant for these parameter values).
>
> Cheers,
> Paul
>
> --
> Paul Abbott Phone: +61 8 6488 2734
> School of Physics, M013 Fax: +61 8 6488 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
>
>