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