Evaluation of the Boys function

$$F_n(x) = \int_0^1 \exp(-xt^2) t^{2n}dt $$

For large $x$, we use the asymptotic form:

$$ F_n(x) \approx \frac{(2n-1)!!}{2^{n+1}} \sqrt{\frac{\pi}{x^{2n+1}}} $$

For small $x$, we Taylor expand about $x = 0$:

$$ F_n(x) \approx \sum_{k=0}^{\infty} \frac{(-x)^k}{k!(2n+2k+1)} $$

For intermediate $x$, we Taylor expand about pretabulated values of $ F_n (x)$:

$$ F_n(x_t+\Delta t) \approx \sum_{k=0}^{\infty} \frac{F_{n+k}(x_t) (-\Delta x)^k}{k!} $$

Boys function for $n=0$:

$$ F_0(x) = erf(\sqrt{x}) \frac{\sqrt{\pi}}{2\sqrt{x}} $$

Boys function for $x=0$:

$$ F_n(0) = \frac{1}{2n+1} $$

Test for $F_0(x)$ for small $x <= 50$:



In [6]:

    
from math import erf
def erfVersionOfBoysFunction(x):
    return erf(sqrt(x))*sqrt(pi)/(2*sqrt(x))



In [24]:

    
s = 0.01; e=50; c = 20
nmax = 100;
xvec = linspace(s,e,c)
F0 = zeros(c)

print "Boys boysF0_small(0);"+ "\n" \
      +"mat F0; rowvec xvec;" + "\n" \
      +"xvec = linspace<rowvec>(" + str(s) + "," + str(e)+ "," + str(c) + ");" \
      + "\n" + "F0 = zeros(" + str(c)+ ",1);"

print "for(uint i = 0; i < " + str(c) +  "; i++){" + "\n" \
      + "  boysF0_small.evaluateBoysFunctions(xvec[i]);" + "\n" \
      + "  F0.row(i) =  boysF0_small.getBoysFunctions();" + "\n" \
      + "}" +"\n"
        
        
for i in range(0, c):
    F0[i]  = erfVersionOfBoysFunction(xvec[i])
    
    print "CHECK_CLOSE(F0(" \
          +str(i)+ ",0)"+ ", " + str("%.16e" %F0[i])\
          + ", 1e-10);"









    



Boys boysF0_small(0);
mat F0; rowvec xvec;
xvec = linspace<rowvec>(0.01,50,20);
F0 = zeros(20,1);
for(uint i = 0; i < 20; i++){
  boysF0_small.evaluateBoysFunctions(xvec[i]);
  F0.row(i) =  boysF0_small.getBoysFunctions();
}

CHECK_CLOSE(F0(0,0), 9.9667664290336333e-01, 1e-10);
CHECK_CLOSE(F0(1,0), 5.3357683580906246e-01, 1e-10);
CHECK_CLOSE(F0(2,0), 3.8551956882369670e-01, 1e-10);
CHECK_CLOSE(F0(3,0), 3.1522027004511449e-01, 1e-10);
CHECK_CLOSE(F0(4,0), 2.7304989839411259e-01, 1e-10);
CHECK_CLOSE(F0(5,0), 2.4424745291117503e-01, 1e-10);
CHECK_CLOSE(F0(6,0), 2.2298057384240719e-01, 1e-10);
CHECK_CLOSE(F0(7,0), 2.0644923632081960e-01, 1e-10);
CHECK_CLOSE(F0(8,0), 1.9312212797082015e-01, 1e-10);
CHECK_CLOSE(F0(9,0), 1.8208209208151274e-01, 1e-10);
CHECK_CLOSE(F0(10,0), 1.7274188563417292e-01, 1e-10);
CHECK_CLOSE(F0(11,0), 1.6470576997839928e-01, 1e-10);
CHECK_CLOSE(F0(12,0), 1.5769603853464209e-01, 1e-10);
CHECK_CLOSE(F0(13,0), 1.5151129820345471e-01, 1e-10);
CHECK_CLOSE(F0(14,0), 1.4600146419605400e-01, 1e-10);
CHECK_CLOSE(F0(15,0), 1.4105209097756166e-01, 1e-10);
CHECK_CLOSE(F0(16,0), 1.3657418098484136e-01, 1e-10);
CHECK_CLOSE(F0(17,0), 1.3249734289737117e-01, 1e-10);
CHECK_CLOSE(F0(18,0), 1.2876507161025572e-01, 1e-10);
CHECK_CLOSE(F0(19,0), 1.2533141373155002e-01, 1e-10);

Test for $F_0(x)$ for large $x > 50$:



In [10]:

    
def factorialBoysFunction(x,n):

    return 1.0/(2**(n+1)) * sqrt(pi/x**(2*n+1))



In [25]:

    
s = 51; e=100; c = 20
nmax = 100;
xvec = linspace(s,e,c)
F0 = zeros(c)

print "Boys boysF0_large(0);"+ "\n" \
      +"xvec = linspace<rowvec>(" + str(s) + "," + str(e)+ "," + str(c) + ");" \
      + "\n" + "F0 = zeros(" + str(c)+ ",1);"

print "for(uint i = 0; i < " + str(c) +  "; i++){" + "\n" \
      + " boysF0_large.evaluateBoysFunctions(xvec[i]);" + "\n" \
      + " F0.row(i) = boysF0_large.getBoysFunctions();" + "\n" \
      + "}" +"\n"
        
        
for i in range(0, c):
    F0[i]  = erfVersionOfBoysFunction(xvec[i])
    
    print "CHECK_CLOSE(F0(" \
          +str(i)+ ",0)"+ ", " + str("%.16e" %F0[i])\
          + ", 1e-10);"









    



Boys boysF0_large(0);
xvec = linspace<rowvec>(51,100,20);
F0 = zeros(20,1);
for(uint i = 0; i < 20; i++){
 boysF0_large.evaluateBoysFunctions(xvec[i]);
 F0.row(i) = boysF0_large.getBoysFunctions();
}

CHECK_CLOSE(F0(0,0), 1.2409659136408727e-01, 1e-10);
CHECK_CLOSE(F0(1,0), 1.2107315290425182e-01, 1e-10);
CHECK_CLOSE(F0(2,0), 1.1826045114934411e-01, 1e-10);
CHECK_CLOSE(F0(3,0), 1.1563509029711123e-01, 1e-10);
CHECK_CLOSE(F0(4,0), 1.1317715652582763e-01, 1e-10);
CHECK_CLOSE(F0(5,0), 1.1086957884293665e-01, 1e-10);
CHECK_CLOSE(F0(6,0), 1.0869762771661817e-01, 1e-10);
CHECK_CLOSE(F0(7,0), 1.0664851770975842e-01, 1e-10);
CHECK_CLOSE(F0(8,0), 1.0471108954515591e-01, 1e-10);
CHECK_CLOSE(F0(9,0), 1.0287555349437555e-01, 1e-10);
CHECK_CLOSE(F0(10,0), 1.0113328058446551e-01, 1e-10);
CHECK_CLOSE(F0(11,0), 9.9476631436519997e-02, 1e-10);
CHECK_CLOSE(F0(12,0), 9.7898814974335960e-02, 1e-10);
CHECK_CLOSE(F0(13,0), 9.6393771031862641e-02, 1e-10);
CHECK_CLOSE(F0(14,0), 9.4956072224460661e-02, 1e-10);
CHECK_CLOSE(F0(15,0), 9.3580841456184838e-02, 1e-10);
CHECK_CLOSE(F0(16,0), 9.2263682201430275e-02, 1e-10);
CHECK_CLOSE(F0(17,0), 9.1000619287057175e-02, 1e-10);
CHECK_CLOSE(F0(18,0), 8.9788048355705335e-02, 1e-10);
CHECK_CLOSE(F0(19,0), 8.8622692545275800e-02, 1e-10);

Testing Boys downward recursion relation:

$$ F_{n-1}(x) = \frac{2xF_n(x)+ \exp(-x)}{2n-1} $$

Test for $F_n(0)$:

Using closed form expression:



In [16]:

    
n_max = 20
F_exact = zeros(n_max+1)

for n in range(0,n_max+1,1):
    F_exact[n] = 1.0/(2*n+1)
    
    print "CHECK_CLOSE(Fn[" + str(n)+ "]"+ ", " +  str("%.16e" % F_exact[n])\
          + ", 1e-10);"









    



CHECK_CLOSE(Fn[0], 1.0000000000000000e+00, 1e-10);
CHECK_CLOSE(Fn[1], 3.3333333333333331e-01, 1e-10);
CHECK_CLOSE(Fn[2], 2.0000000000000001e-01, 1e-10);
CHECK_CLOSE(Fn[3], 1.4285714285714285e-01, 1e-10);
CHECK_CLOSE(Fn[4], 1.1111111111111110e-01, 1e-10);
CHECK_CLOSE(Fn[5], 9.0909090909090912e-02, 1e-10);
CHECK_CLOSE(Fn[6], 7.6923076923076927e-02, 1e-10);
CHECK_CLOSE(Fn[7], 6.6666666666666666e-02, 1e-10);
CHECK_CLOSE(Fn[8], 5.8823529411764705e-02, 1e-10);
CHECK_CLOSE(Fn[9], 5.2631578947368418e-02, 1e-10);
CHECK_CLOSE(Fn[10], 4.7619047619047616e-02, 1e-10);
CHECK_CLOSE(Fn[11], 4.3478260869565216e-02, 1e-10);
CHECK_CLOSE(Fn[12], 4.0000000000000001e-02, 1e-10);
CHECK_CLOSE(Fn[13], 3.7037037037037035e-02, 1e-10);
CHECK_CLOSE(Fn[14], 3.4482758620689655e-02, 1e-10);
CHECK_CLOSE(Fn[15], 3.2258064516129031e-02, 1e-10);
CHECK_CLOSE(Fn[16], 3.0303030303030304e-02, 1e-10);
CHECK_CLOSE(Fn[17], 2.8571428571428571e-02, 1e-10);
CHECK_CLOSE(Fn[18], 2.7027027027027029e-02, 1e-10);
CHECK_CLOSE(Fn[19], 2.5641025641025640e-02, 1e-10);
CHECK_CLOSE(Fn[20], 2.4390243902439025e-02, 1e-10);

Test for $F_n(x)$ for small $x <= 50$ , when $n \rightarrow 0$:

Using error function:



In [28]:

    
s = 0.01; e=50; c = 20
xvec = linspace(s,e,c)
F0 = zeros(c)
n = 15
print "Boys boysF0_small("+str(n)+");"+ "\n" \
      +"mat F0; rowvec xvec;" + "\n" \
      +"xvec = linspace<rowvec>(" + str(s) + "," + str(e)+ "," + str(c) + ");" \
      + "\n" + "F0 = zeros(" + str(c)+ "," + str(n+1)+ ");"

print "for(uint i = 0; i < " + str(c) +  "; i++){" + "\n" \
      + "  boysF0_small.evaluateBoysFunctions(xvec[i]);" + "\n" \
      + "  F0.row(i) =  boysF0_small.getBoysFunctions();" + "\n" \
      + "}" +"\n"
        
        
for i in range(0, c):
    F0[i]  = erfVersionOfBoysFunction(xvec[i])
    
    print "CHECK_CLOSE(F0(" \
          +str(i)+ ",0)"+ ", " + str("%.16e" % F0[i])\
          + ", 1e-10);"









    



Boys boysF0_small(15);
mat F0; rowvec xvec;
xvec = linspace<rowvec>(0.01,50,20);
F0 = zeros(20,16);
for(uint i = 0; i < 20; i++){
  boysF0_small.evaluateBoysFunctions(xvec[i]);
  F0.row(i) =  boysF0_small.getBoysFunctions();
}

CHECK_CLOSE(F0(0,0), 9.9667664290336333e-01, 1e-10);
CHECK_CLOSE(F0(1,0), 5.3357683580906246e-01, 1e-10);
CHECK_CLOSE(F0(2,0), 3.8551956882369670e-01, 1e-10);
CHECK_CLOSE(F0(3,0), 3.1522027004511449e-01, 1e-10);
CHECK_CLOSE(F0(4,0), 2.7304989839411259e-01, 1e-10);
CHECK_CLOSE(F0(5,0), 2.4424745291117503e-01, 1e-10);
CHECK_CLOSE(F0(6,0), 2.2298057384240719e-01, 1e-10);
CHECK_CLOSE(F0(7,0), 2.0644923632081960e-01, 1e-10);
CHECK_CLOSE(F0(8,0), 1.9312212797082015e-01, 1e-10);
CHECK_CLOSE(F0(9,0), 1.8208209208151274e-01, 1e-10);
CHECK_CLOSE(F0(10,0), 1.7274188563417292e-01, 1e-10);
CHECK_CLOSE(F0(11,0), 1.6470576997839928e-01, 1e-10);
CHECK_CLOSE(F0(12,0), 1.5769603853464209e-01, 1e-10);
CHECK_CLOSE(F0(13,0), 1.5151129820345471e-01, 1e-10);
CHECK_CLOSE(F0(14,0), 1.4600146419605400e-01, 1e-10);
CHECK_CLOSE(F0(15,0), 1.4105209097756166e-01, 1e-10);
CHECK_CLOSE(F0(16,0), 1.3657418098484136e-01, 1e-10);
CHECK_CLOSE(F0(17,0), 1.3249734289737117e-01, 1e-10);
CHECK_CLOSE(F0(18,0), 1.2876507161025572e-01, 1e-10);
CHECK_CLOSE(F0(19,0), 1.2533141373155002e-01, 1e-10);

Test for $F_n(x)$ for large $x > 50$, when $n \rightarrow 0$:

Using error function:



In [29]:

    
s = 51; e=100; c = 20
xvec = linspace(s,e,c)
F0 = zeros(c)
n = 15
print "Boys boysF0_large("+str(n)+");"+ "\n" \
      +"xvec = linspace<rowvec>(" + str(s) + "," + str(e)+ "," + str(c) + ");" \
      + "\n" + "F0 = zeros(" + str(c)+ ","+str(n+1)+");"

print "for(uint i = 0; i < " + str(c) +  "; i++){" + "\n" \
      + " boysF0_large.evaluateBoysFunctions(xvec[i]);" + "\n" \
      + " F0.row(i) = boysF0_large.getBoysFunctions();" + "\n" \
      + "}" +"\n"
        
        
for i in range(0, c):
    F0[i]  = erfVersionOfBoysFunction(xvec[i])
    
    print "CHECK_CLOSE(F0(" \
          +str(i)+ ",0)"+ ", " + str("%.16e" %F0[i])\
          + ", 1e-10);"









    



Boys boysF0_large(15);
xvec = linspace<rowvec>(51,100,20);
F0 = zeros(20,16);
for(uint i = 0; i < 20; i++){
 boysF0_large.evaluateBoysFunctions(xvec[i]);
 F0.row(i) = boysF0_large.getBoysFunctions();
}

CHECK_CLOSE(F0(0,0), 1.2409659136408727e-01, 1e-10);
CHECK_CLOSE(F0(1,0), 1.2107315290425182e-01, 1e-10);
CHECK_CLOSE(F0(2,0), 1.1826045114934411e-01, 1e-10);
CHECK_CLOSE(F0(3,0), 1.1563509029711123e-01, 1e-10);
CHECK_CLOSE(F0(4,0), 1.1317715652582763e-01, 1e-10);
CHECK_CLOSE(F0(5,0), 1.1086957884293665e-01, 1e-10);
CHECK_CLOSE(F0(6,0), 1.0869762771661817e-01, 1e-10);
CHECK_CLOSE(F0(7,0), 1.0664851770975842e-01, 1e-10);
CHECK_CLOSE(F0(8,0), 1.0471108954515591e-01, 1e-10);
CHECK_CLOSE(F0(9,0), 1.0287555349437555e-01, 1e-10);
CHECK_CLOSE(F0(10,0), 1.0113328058446551e-01, 1e-10);
CHECK_CLOSE(F0(11,0), 9.9476631436519997e-02, 1e-10);
CHECK_CLOSE(F0(12,0), 9.7898814974335960e-02, 1e-10);
CHECK_CLOSE(F0(13,0), 9.6393771031862641e-02, 1e-10);
CHECK_CLOSE(F0(14,0), 9.4956072224460661e-02, 1e-10);
CHECK_CLOSE(F0(15,0), 9.3580841456184838e-02, 1e-10);
CHECK_CLOSE(F0(16,0), 9.2263682201430275e-02, 1e-10);
CHECK_CLOSE(F0(17,0), 9.1000619287057175e-02, 1e-10);
CHECK_CLOSE(F0(18,0), 8.9788048355705335e-02, 1e-10);
CHECK_CLOSE(F0(19,0), 8.8622692545275800e-02, 1e-10);



In [8]:



In [8]:



In [8]:



In [8]:



In [8]:

Test for $F_0(x)$ for small $x < 0.09$:



In [9]:

    
s = 0.01; e=0.09; c = 10
nmax = 100;
xvec = linspace(s,e,c)
Fpvec = zeros(c)
F_exact = zeros(c)

print "vec xvec, Fpvec;" + "\n" \
      +"xvec = linspace(" + str(s) + "," + str(e)+ "," + str(c) + ");" \
      + "\n" + "Fpvec = zeros(" + str(c)+ ");"

print "for(uint i = 0; i < " + str(c) +  "; i++){" + "\n" \
        "  for(uint n = " + str(nmax)+ "; n > 0; n--){" + "\n" \
      + "           Fpvec[i] = boys.downwardRecursion(xvec[i],n);" + "\n" \
      + "}" +"\n" \
      + "  }"
        
        
for i in range(0, c):
    F_exact[i]  = erfVersionOfBoysFunction(xvec[i])
    Fpvec[i] = factorialBoysFunction(xvec[i], 0)
    
    print "CHECK_CLOSE(Fpvec[" \
          +str(i)+ "]"+ ", " + str(F_exact[i])\
          + ", 1e-10);"









    



vec xvec, Fpvec;
xvec = linspace(0.01,0.09,10);
Fpvec = zeros(10);
for(uint i = 0; i < 10; i++){
  for(uint n = 100; n > 0; n--){
           Fpvec[i] = boys.downwardRecursion(xvec[i],n);
}
  }
CHECK_CLOSE(Fpvec[0], 0.996676642903, 1e-10);
CHECK_CLOSE(Fpvec[1], 0.993739222842, 1e-10);
CHECK_CLOSE(Fpvec[2], 0.990817393658, 1e-10);
CHECK_CLOSE(Fpvec[3], 0.987911056819, 1e-10);
CHECK_CLOSE(Fpvec[4], 0.985020114475, 1e-10);
CHECK_CLOSE(Fpvec[5], 0.982144469445, 1e-10);
CHECK_CLOSE(Fpvec[6], 0.97928402522, 1e-10);
CHECK_CLOSE(Fpvec[7], 0.976438685953, 1e-10);
CHECK_CLOSE(Fpvec[8], 0.973608356454, 1e-10);
CHECK_CLOSE(Fpvec[9], 0.97079294219, 1e-10);

Test for $F_0(x)$ for large $x > 20$:



In [10]:

    
s = 20; e=100; c = 10
nmax = 10;
xvec = linspace(s,e,c)
Fpvec = zeros(c)
F_exact = zeros(c)

print "vec xvec, Fpvec;" + "\n" \
      +"xvec = linspace(" + str(s) + "," + str(e)+ "," + str(c) + ");" \
      + "\n" + "Fpvec = zeros(" + str(c)+ ");"

print "for(uint i = 0; i < " + str(c) +  "; i++){" + "\n" \
        "  for(uint n = " + str(nmax)+ "; n > 0; n--){" + "\n" \
      + "           Fpvec[i] = boys.downwardRecursion(xvec[i],n);" + "\n" \
      + "}" +"\n" \
      + "  }"
        
        
for i in range(0, c):
    F_exact[i]  = erfVersionOfBoysFunction(xvec[i])
    Fpvec[i] = factorialBoysFunction(xvec[i], 0)
    
    print "CHECK_CLOSE(Fpvec[" \
          +str(i)+ "]"+ ", " + str(F_exact[i])\
          + ", 1e-10);"









    



vec xvec, Fpvec;
xvec = linspace(20,100,10);
Fpvec = zeros(10);
for(uint i = 0; i < 10; i++){
  for(uint n = 10; n > 0; n--){
           Fpvec[i] = boys.downwardRecursion(xvec[i],n);
}
  }
CHECK_CLOSE(Fpvec[0], 0.19816636483, 1e-10);
CHECK_CLOSE(Fpvec[1], 0.164884382227, 1e-10);
CHECK_CLOSE(Fpvec[2], 0.144187209502, 1e-10);
CHECK_CLOSE(Fpvec[3], 0.12973033818, 1e-10);
CHECK_CLOSE(Fpvec[4], 0.118899818928, 1e-10);
CHECK_CLOSE(Fpvec[5], 0.110395710425, 1e-10);
CHECK_CLOSE(Fpvec[6], 0.103489008863, 1e-10);
CHECK_CLOSE(Fpvec[7], 0.0977350491129, 1e-10);
CHECK_CLOSE(Fpvec[8], 0.0928451600494, 1e-10);
CHECK_CLOSE(Fpvec[9], 0.0886226925453, 1e-10);



In [10]: