我正在阅读“Fortran 77 中的数值方法：科学计算的艺术”（第二版）一书，我遇到了一些一维函数数值积分的方法。更具体地说，Gauss-Laguerre、Gauss Hermite 和 Gauss Jacobi 权重和横坐标吸引了我。我需要它们，因为我有一些行为不端的被积函数，我需要对它们进行数值积分。据我所知，这些例程没有python实现，最接近这些方法的是http://docs.scipy.org/doc/numpy/reference/generated/numpy.polynomial.hermite.hermgauss.html

但我觉得文档有点不完整。我不确定如何使用这种方法，并且有点困惑为什么在 python 中这些类型的求积没有标准方法。

编辑：

我已经在示例问题上实现并比较了 scipy 正交与 Gauss-Hermite 正交：

\begin{aligned} I = \int_{- \infty}^{\infty} (i ζ (t) + \bar{i ζ (t)}) \exp (i ω t) d t \end{aligned}

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

这里 $i=\sqrt{-1}$ 和 $\zeta(t)$ 是一个复杂的函数 $t$ . overbear 表示复共轭。 $\zeta(t)$ 不是标准函数，它是我计算出来的，样本可以在http://speedy.sh/SX6JU/seData.npy作为 numpy 数组找到。数组中的每一行的形式为 [t real(zeta) imag(zeta)]，即函数在时间 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)

在此处输入图像描述

我不确定正确的结果应该是什么，但我相信结果中不应该存在“噪音”。

相反，我尝试了一个 Gauss-Hermite 例程：

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)

在此处输入图像描述

我相信这更符合我应该得到的，但是我不确定我做对了。似乎 Gauss-hermite 的阶数非常高，而在低阶（90）看起来我的结果显示了某种类型的周期性/混叠。

Python 中是否有 Gauss-Laguerre 积分例程？

编辑：