Re: Numerical accuracy of Hypergeometric2F1
- To: mathgroup at smc.vnet.net
- Subject: [mg55857] Re: Numerical accuracy of Hypergeometric2F1
- From: Paul Abbott <paul at physics.uwa.edu.au>
- Date: Fri, 8 Apr 2005 01:36:29 -0400 (EDT)
- Organization: The University of Western Australia
- References: <d3043c$nei$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
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