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:
It's a combination of:
And:
And:
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.