MathGroup Archive 2005

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

Search the Archive

Re: Numerical accuracy of Hypergeometric2F1

  • To: mathgroup at smc.vnet.net
  • Subject: [mg55803] Re: Numerical accuracy of Hypergeometric2F1
  • From: "Christos Argyropoulos M.D." <chrisarg at fuse.net>
  • Date: Wed, 6 Apr 2005 03:12:02 -0400 (EDT)
  • Sender: owner-wri-mathgroup at wolfram.com

Hello,
Thanx everyone for pointing out the subtleties of making Mathematica 
understand precision of input (I always thought that SetPrecision threads 
over its first input). 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).
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]) 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

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

Cheers and thanks everyone for sharing how to get the best out of the 
arbirtrary-precision point arithmetic with Mathematica

Christos Argyropoulos MD PhD


Mathematica Code
--------------------------------------------------------------------------------------------------------------------
In[1]:=
\!\(fpol[x_, k_] :=
    Module[{}, coeff = Table[1, {k + 1}];
      For[i = 2, i \[LessEqual] k + 1, \(i++\),
        coeff[\([i]\)] =
          coeff[\([i -
                  1]\)]\ \(\(-k\) + i - 1\)\/\(i - k - 1/2\)]; \
{Hypergeometric2F1[k + 1/2, 1, 3/2,
          x], \(1\/\(2  k -
                1\)\) \(\[Sum]\+\(i = 1\)\%\(k + 1\)\((\((1 - x)\)\^\(-i\)\ 
\
coeff[\([i]\)])\)\)}]\[IndentingNewLine]
  fhorn[x_, k_] :=
    Module[{}, coeff = Table[1, {k + 1}];
      For[i = 2, i \[LessEqual] k + 1, \(i++\),
        coeff[\([i]\)] =
          coeff[\([i -
                  1]\)]\ \(\(-k\) + i - 1\)\/\(i - k - 1/2\)]; \
{Hypergeometric2F1[k + 1/2, 1, 3/2, x], \(1\/\(2  k - 1\)\)
          Module[{y = 1/\((1 - x)\), res = 1}, res = 0;
            For[i = k + 1, i \[GreaterEqual] \ 1\ , \(i--\),
              res = \((res + coeff[\([i]\)])\)\ y\ ]; res]}]\)
In[3]:=
x=9/10;
k=90.;
workprec=16;
printprec=16;
In[7]:=

Print["Workprec = ",workprec];
fpol[x,k]//SetPrecision[#,printprec]&
fpol[N[x],k]//SetPrecision[#,16]&
fpol[SetPrecision[x,workprec],k]//SetPrecision[#,printprec]&
fpol[SetPrecision[x,workprec],SetPrecision[k,workprec]]//
  SetPrecision[#,printprec]&
fpol[x,Rationalize[k]]//SetPrecision[#,printprec]&
workprec=20;
Print["Workprec = ",workprec];
fpol[x,k]//SetPrecision[#,printprec]&
fpol[N[x],k]//SetPrecision[#,16]&
fpol[SetPrecision[x,workprec],k]//SetPrecision[#,printprec]&
fpol[SetPrecision[x,workprec],SetPrecision[k,workprec]]//
  SetPrecision[#,printprec]&
fpol[x,Rationalize[k]]//SetPrecision[#,printprec]&
workprec=25;
Print["Workprec = ",workprec];
fpol[x,k]//SetPrecision[#,printprec]&
fpol[N[x],k]//SetPrecision[#,16]&
fpol[SetPrecision[x,workprec],k]//SetPrecision[#,printprec]&
fpol[SetPrecision[x,workprec],SetPrecision[k,workprec]]//
  SetPrecision[#,printprec]&
fpol[x,Rationalize[k]]//SetPrecision[#,printprec]&
From In[7]:=
Workprec = \[InvisibleSpace]16
Out[8]=
\!\({7.20338647551431814577488138862587`15.9546*^95,
    9.860651611217270453862251879131`15.9546*^88}\)
Out[9]=
\!\({7.20338647551431814577488138862587`15.9546*^95,
    9.860651611217469754211158880883`15.9546*^88}\)
Out[10]=
\!\({7.20338647551431814577488138862587`15.9546*^95,
    9.860651611217466927255855235469`15.9546*^88}\)
Out[11]=
\!\({9.990653929442517915012215795621`15.9546*^88,
    9.860651611217461201585792023911`15.9546*^88}\)
Out[12]=
\!\({9.860651611217264799951644588301`15.9546*^88,
    9.860651611217264799951644588301`15.9546*^88}\)
From In[7]:=
Workprec = \[InvisibleSpace]20
Out[15]=
\!\({7.20338647551431814577488138862587`15.9546*^95,
    9.860651611217270453862251879131`15.9546*^88}\)
Out[16]=
\!\({7.20338647551431814577488138862587`15.9546*^95,
    9.860651611217469754211158880883`15.9546*^88}\)
Out[17]=
\!\({7.20338647551431814577488138862587`15.9546*^95,
    9.860651611217270453862251879131`15.9546*^88}\)
Out[18]=
\!\({9.8606020395819286190853862`15.9546*^88,
    9.860651611217264267879523610443`15.9546*^88}\)
Out[19]=
\!\({9.860651611217264799951644588301`15.9546*^88,
    9.860651611217264799951644588301`15.9546*^88}\)
From In[7]:=
Workprec = \[InvisibleSpace]25
Out[22]=
\!\({7.20338647551431814577488138862587`15.9546*^95,
    9.860651611217270453862251879131`15.9546*^88}\)
Out[23]=
\!\({7.20338647551431814577488138862587`15.9546*^95,
    9.860651611217469754211158880883`15.9546*^88}\)
Out[24]=
\!\({7.20338647551431814577488138862587`15.9546*^95,
    9.860651611217270453862251879131`15.9546*^88}\)
Out[25]=
\!\({9.8606516120571455699381674`15.9546*^88,
    9.8606516112172642678795376`15.9546*^88}\)
Out[26]=
\!\({9.860651611217264799951644588301`15.9546*^88,
    9.860651611217264799951644588301`15.9546*^88}\)
In[27]:=
workprec=25;
x=9/10;
k=10.;
Print["Workprec = ",workprec];
fpol[x,k]//SetPrecision[#,printprec]&
fpol[N[x],k]//SetPrecision[#,16]&
fpol[SetPrecision[x,workprec],k]//SetPrecision[#,printprec]&
fpol[SetPrecision[x,workprec],SetPrecision[k,workprec]]//
  SetPrecision[#,printprec]&
fpol[x,Rationalize[k]]//SetPrecision[#,printprec]&
From In[27]:=
Workprec = \[InvisibleSpace]25
Out[31]=
\!\({2.991232093254027843475341796875`15.9546*^9,
    2.991232093254022121429443359375`15.9546*^9}\)
Out[32]=
\!\({2.991232093254027843475341796875`15.9546*^9,
    2.9912320932540283203125`15.9546*^9}\)
Out[33]=
\!\({2.991232093254027843475341796875`15.9546*^9,
    2.991232093254022121429443359375`15.9546*^9}\)
Out[34]=
\!\({2.9912320932540215202753902439`15.9546*^9,
    2.9912320932540215202753902444`15.9546*^9}\)
Out[35]=
\!\({2.99123209325402164459228515625`15.9546*^9,
    2.99123209325402164459228515625`15.9546*^9}\)
------------------------------------------------------------------------------------------------

Fortran Implementation
source for the hypergeometric function : 
http://cfa-www.harvard.edu/~rcf/codes/hyp.f
test file source code:

c********************************************************************
c
c   program :  hyptest

c***********************************************************************
c

      real*8 a,b,c,z,w,re,im
      a=1.d0
      b=90.5d0
      c=1.5d0
      z=0.9d0
      call hyp(z,a,b,c,re,im)
 write(*,900)re
900   format (g22.16)
1000  stop
      end

compile as:
   g77  hyp.f  hyptest.f 



  • Prev by Date: Re: Re: Mathematica bug in handling trigonometric functions? (and more)
  • Next by Date: ploting functions slowly
  • Previous by thread: Re: Re: Re[2]: Re: Numerical accuracy of Hypergeometric2F1
  • Next by thread: Re: Numerical accuracy of Hypergeometric2F1