Is there a Gauss-Laguerre integration routine in Python?

I am reading the book "Numerical Recipes in Fortran 77: The Art of Scientific Computing" (Second Edition) and I came across some methods for numerical integration of 1D functions. More specifically the Gauss-Laguerre, Gauss Hermite, and Gauss Jacobi weights and abscissas appealed to me. I need them because I have some ill-behaved integrands that I need to integrate numerically. As far as I can tell there is no python implementation of these routines, the closest to these methods was this http://docs.scipy.org/doc/numpy/reference/generated/numpy.polynomial.hermite.hermgauss.html

but I feel that documentation is a bit incomplete. I'm not sure what to do with this method and slightly puzzled why is there no standard method for these types of quadrature in python.

EDIT:

I have implemented and compared scipy quadrature with Gauss-Hermite Quadrature on the example problem:

\begin{align} I = \int_{-\infty}^{\infty} \left( i\zeta(t) + \overline{i\zeta(t)} \right) \exp(i\omega t) \, \mathrm{d}t \end{align}

here $i=\sqrt{-1}$ and $\zeta(t)$ is a complex function of $t$. The overbear denotes complex conjugate. $\zeta(t)$ is not a standard function, it is something I calculated and the samples can be found at http://speedy.sh/SX6JU/seData.npy as a numpy array. Each row in the array is of form [t real(zeta) imag(zeta)] i.e. the sampled real and imaginary values of the function at time t.

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate

def integrand(t, omega):
x1 = np.interp(t, T1[:,0], T1[:,1])
x2 = np.interp(t, T1[:,0], T1[:,2])
z0 = x1+x2*1j
zeta = z0/np.sqrt(z0**2+1)
return (1j*zeta+np.conjugate(1j*zeta))*np.exp(1j*omega*t)

def integral(omega):
def real_func(x,omega):
return np.real(integrand(x,omega))
def imag_func(x,t):
return np.imag(integrand(x,omega))
real_integral = integrate.quad(real_func, -80, 80, args=(omega))
imag_integral = integrate.quad(imag_func, -80, 80, args=(omega))
return real_integral[0] + 1j*imag_integral[0]

vintegral = np.vectorize(integral)

omega = np.linspace(-30, 30, 601)
I = vintegral(omega)


Im not sure what the correct result should be but the "noisiness" present in the result should not be there I believe.

Instead I tried a Gauss-Hermite routine:

ti_her, wi_her, = np.polynomial.hermite.hermgauss(90)

vintegrand = np.vectorize(integrand)

def integral_hermite(omega):
real_integral = np.dot(vintegrand(ti_her, omega).real, wi_her)
imag_integral = np.dot(vintegrand(ti_her, omega).imag, wi_her)
return real_integral + 1j*imag_integral

vintegral_hermite = np.vectorize(integral_hermite)

I2 = vintegral_hermite(omega)


I believe this is more in line with what I should be getting , however I am not sure I have done the correct thing. It seems the order of Gauss-hermite is very high and at low orders (90) it looks like my result shows some type of periodicity/aliasing.

There are standard methods for these types of quadrature in Python, in NumPy and SciPy:

and other routines.

Given a degree of quadrature (essentially, the degree of the interpolant used to approximate the integrand), these Gauss-type quadrature functions return quadrature points $x_{i}$ and weight functions $w_{i}$ so that you can approximate the integral

\begin{align} \int_{a}^{b}f(x)\,\mathrm{d}x \end{align}

for appropriate $a$ and $b$ by

\begin{align} \sum_{i=0}^{N}f(x_{i})w_{i}; \end{align}

it's up to you to use the quadrature points and weights in a summation to calculate an approximation to your integral.

• Thanks for your answer. I get a problem with the Gauss-Laguerre routine in that the some of the $x_i$ returned for a 175 degree are above the range for which I have sampled values of $f$. Does that mean I need to compute $f$ at more points in order to obtain 175 point accuracy? I just seems strange that only a single integer determines the sample points and weights and does not take into consideration the behaviour of the integrand. Sep 21, 2014 at 10:56
• @Jack, degree 175 is too high. Way too high. This is equivalent to fitting a 175th degree polynomial through your function and integrating that. $x^{175} + c$ is a very hard thing for a computer to compute with any level of accuracy in finite-precision arithmetic. You will need to divide your interval of integration in to man subintervals and integrate over those with a lower degree quadrature rule. One of the adaptive routines that Geoff linked to is probably the best way to do this, but you can also do it by hand. Sep 21, 2014 at 13:59
• @BillBarth ya now that you mention it I understand that causes some problems. But I'm following a paper which claims to have used a 175 point Gauss Laguerre routine. Could they have been referring to something else? Sep 21, 2014 at 15:42
• @Jack, probably not. Maybe you could update the post with the paper you are trying to follow. Sep 21, 2014 at 16:07
• @BillBarth Please see my edit. I have now implemented and compared a Gauss-Hermite routine with 175 points. I would very much appreciate if you could have a look at it. As for the paper, I only have it on hard copy. I can scan relevant parts of it but it does say what I previously stated. Thanks! Sep 22, 2014 at 14:56