4
$\begingroup$

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

T1 =  np.load('seData.npy') #: [Time, x1, x2] from downloaded file


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)

enter image description here

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)

enter image description here

enter image description here

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.

$\endgroup$
6
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – Dipole Sep 21 '14 at 10:56
  • 2
    $\begingroup$ @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. $\endgroup$ – Bill Barth Sep 21 '14 at 13:59
  • $\begingroup$ @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? $\endgroup$ – Dipole Sep 21 '14 at 15:42
  • $\begingroup$ @Jack, probably not. Maybe you could update the post with the paper you are trying to follow. $\endgroup$ – Bill Barth Sep 21 '14 at 16:07
  • $\begingroup$ @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! $\endgroup$ – Dipole Sep 22 '14 at 14:56
2
$\begingroup$

quadpy (a small project of mine) might be helpful. For quadrature of line segments, it features

  • Chebyshev-Gauß,
  • Clenshaw-Curtis,
  • Fejér (both types),
  • Gauß-Hermite,
  • Gauß-Laguerre,
  • Gauß-Legendre,
  • Gauß-Lobatto,
  • Gauß-Patterson (7 schemes up to degree 191),
  • Gauß-Radau,
  • open and closed Newton-Cotes.
$\endgroup$

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.