# 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 + 1j*imag_integral

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. – Dipole Sep 21 '14 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. – Bill Barth Sep 21 '14 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? – Dipole Sep 21 '14 at 15:42
• @Jack, probably not. Maybe you could update the post with the paper you are trying to follow. – Bill Barth Sep 21 '14 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! – Dipole Sep 22 '14 at 14:56