Welcome, guest | Sign In | My Account | Store | Cart

A new method for computing the hypergeometric function 1F1(a,b,t) (Python recipe) by Fernando Nieuwveldt
ActiveState Code (http://code.activestate.com/recipes/577133/)

I present a method of computing the 1F1(a,b,x) function using a contour integral. The method is based on a numerical inversion, basically the Laplace inversion. Integral is 1F1(a,b,x) = Gamma(b)/2\pi i \int_\rho exp(zx)z^(-b)(1+x/z)^(-a)dz, \rho is taken as a Talbot contour. The Talbot method is applied with the use of the midpoint rule for numerical integration. Here the user must give the number of function evaluations and this may vary from problem to problem. It is very easy to implement with only a few lines of code and it is very accurate even for large arguments.

Python, 54 lines
Download
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
################################################################################
## Created by Fernando Damian Nieuwveldt                                        
## Date : 26 October 2009                                                      
## email : fdnieuwveldt@gmail.com
## This was part work of my masters thesis in Applied Mathematics at the
## University of Stellenbosch, South Africa
## Thesis title : A Survey of Computational Methods for Pricing Asian Options
## For reference details contact me via email.
################################################################################
# Example : 
# >>>Hypergeometric1F1(5,2,100-1000j,50)
# (7.0028644420385242e+050+8.9737757674587575e+050j)
#
# using mpmath
# Hypergeometric1F1(-1000,1,1000,2000)
# mpc(real='-2.5938207833620713892854e+215', imag='-1.1411010130233474563877e+202')
#
# Other references :
# A. Erdélyi, W. Magnus, F. Oberhettinger, and F. Tricomi. Higher
# transcendental functions. Vol. I. Robert E. Krieger Publishing Co.
# Inc., Melbourne, Fla., 1981, pg 273
#
# L.N.Trefethen, J.A.C.Weideman, and T.Schmelzer. Talbot quadratures
# and rational approximations. BIT. Numerical Mathematics,
# 46(3):653 670, 2006.


from mpmath import mp,mpf,mpc,pi,sin,tan,exp,gamma
mp.dps = 100

def F(a,b,t,s):
   
    return s**(-b)*(1+1.0/s)**(a-b)

def Hypergeometric1F1(a,b,t,N):
    
#   Initiate the stepsize
    h = 2*pi/N;  
    
    c1 = mpf('0.5017')
    c2 = mpf('0.6407')
    c3 = mpf('0.6122')
    c4 = mpc('0','0.2645')
    
#   The for loop is evaluating the Laplace inversion at each point theta
#   which is based on the trapezoidal rule
    ans = 0.0
    for k in range(N):
      theta = -pi + (k+0.5)*h
      z    = N/t*(c1*theta/tan(c2*theta) - c3 + c4*theta)
      dz   = N/t*(-c1*c2*theta/sin(c2*theta)**2 + c1/tan(c2*theta)+c4)
      ans += exp(z*t)*F(a,b,t,z)*dz
          
    return gamma(b)*t**(1-b)*exp(t)*((h/(2j*pi))*ans)
Tags: function, hypergeometric, laplace_inversion

3 comments

Dieter Kadelka 14 years, 1 month ago  # | flag

Hello, The code is interesting, but the results are not correct. Omit line 38 and try outside your program mp.dps = 500. Then the result of Hypergeometric1F1(-1000,1,1000,2000) is (after some time) mpc(real='-259382078336200571793976408157920288042214518531712673813454643437412674426207965313223390138295462654426024807574715475989846138108832184793462174037341015520450763218289713126187040486437468492716746181595162186973.12131666130303488229807541771487147053002114809544835515314563249372531758751614602957353905062735276883204624477834916510930414363243787643352987882001313314012388126921567916365588372887364753590544416460014775157287991471908599653038787056232234618336392000611051724230417748590597499', imag='-1.5718096942928172818623515911296029575960431665555802398073479457533119211090161906814926041334408542490957098973895163890077502573645139475267240175982017648898833360300439873436488063125560789189919473997102862810972899232104169686340314592307414858060039198887212982892908387627919450508644705336452713536114884280606969386659775902646727628160796097999063335649324849133354670290540738132814392818146019760591452841018302838015376114300145875425801137921458861778563896774083610331655099595828389592e-278')

which seems to be almost the correct value (computed with mpmath-0.14) hyp1f1(-1000,1,1000) = -2.5938207833620057179397640815792028804221451853171285e+215

What is missing is an automatic precision control.

Fernando Nieuwveldt (author) 14 years, 1 month ago  # | flag

Hi Dieter. I made the correction you suggested. Using mp.dps = 100 I got the following error

>>> error=(hyp1f1(-1000,1,1000)-Hypergeometric1F1(-1000,1,1000,2000))/hyp1f1(-1000,1,1000)
>>> error
-4.26638919839693558720148557862920528979634425017803934138236098443677542785701292944517572090479952751e-94

Also when the values a,b,t are real one can use symmetry as 

>>> for k in range(N/2):

and then need to multiply the final answer with 2:

>>> return 2*gamma(b)*t**(1-b)*exp(t)*((h/(2j*pi))*ans)

For a fast and accurate method using floating point arithmetic this integral approach looks like a good one, but the need of a truncation error comes to light.

Dieter Kadelka 14 years, 1 month ago  # | flag

Hi Fernando, there are still great problems with your code otherwise. Try Hypergeometric1F1(-1000,1,1000,5000). The result of this computation is useless. N = 2000 gives the correct result (N = 1400 seems to be good also), but for larger N the rounding errors are becoming larger and larger. Without automatic precision control and when the result is not known: What N should be chosen and what mp.dps? mp.dps depends on N, that's the problem. Convergence with mp.dps fixed can not be guaranteed.