Gauss-Legendre over intervals -x -> infinity: adaptive algorithm to transform weights and nodes efficiently

I think that code does the job:

import numpy as np 
import math

deg = 10
x, w = np.polynomial.legendre.leggauss(deg)

def function(x): 
    # the function to integrate
    return math.exp(-x)

def function2(x, a): 
    return function(a+x/(1-x))/((1-x)**2); 

def anotherOne(x, a): 
    return 0.5 * function2(x/2 + 1/2, a)

def integrate(deg, a): 
    sum = 0 
    x, w = np.polynomial.legendre.leggauss(deg)
    for i in range(deg): 
        print("sum({}) += {} * {} (eval in {})".format(sum, w[i], anotherOne(x[i], a), x[i]))
        sum += w[i]*anotherOne(x[i], a)
    return sum; 

print("result"); 
print(integrate(10, 1))

It combines your equation to integrate from a to inf and the equation to change the bounds of an integral.

I hope it solves your problem (it works for exp(-x) at least) :)

If you want an inline computation, the program does the sum of: enter image description here

It's a combination of:

enter image description here

And:

enter image description here

And:

enter image description here


In "Numerical Programming: A Practical Guide for Scientists and Engineers Using Python and C/C++" by Titus A. Beu, you can find the methods in the code samples integral.py and specfunc.py here: http://phys.ubbcluj.ro/~tbeu/INP/libraries.html You call the function xGaussLag(a, deg) which calls Laguerre from the other .py file and returns your adjusted (x,w) between a and infinity. Here's how to set this up (note just above deg=80 it is very slow, I'm just showing you how to apply it by modifying lines above):

x, w = np.array(xGaussLag(a,deg))
gauss = sum(w * integrand(x, flag, F, K, vol, T2, T1))

Gets pretty close convergence on deg=80 (faster) but I just put the eps=1e-13 in xGaussLag and pushed the deg=150 with these results, nonetheless faster than quad by 33%:

The QUADPACK solution: 0.149221620346 with error: 1.49870924498e-12 Gauss-Legendre solution: 0.149238273747 Difference between QUADPACK and Gauss-Legendre: 1.66534003601e-05

In Cython this is 6x faster than straight Python BTW still too slow so I'm going to try the "FastGL" package with the answer from @Alexis for now, just posting as I think this will be useful for other SO users in the future.