python之Scipy:Hermite 函数与正交权重的集成

我想使用 scipy.integrate.quad 对两个时移和频移的 Hermite 函数的乘积进行积分。

但是，由于包含大阶多项式，因此会出现数值错误。这是我的代码:

import numpy as np 
import scipy.integrate 
import scipy.special as sp 
from math import pi 
 
 
def makeFuncs(): 
    # Create the 0th, 4th, 8th, 12th and 16th order hermite function 
    return [lambda t, n=n: np.exp(-0.5*t**2)*sp.hermite(n)(t) for n in np.arange(5)*4] 
 
def ambgfun(funcs, i, k, tau, f): 
    # Integrate f1(t)*f2(t+tau)*exp(-j2pift) over t from -inf to inf 
    f1 = funcs[i] 
    f2 = funcs[k] 
    func = lambda t: np.real(f1(t) * f2(t+tau) * np.exp(-1j*(2*pi)*f*t)) 
    return scipy.integrate.quad(func, -np.inf, np.inf) 
 
def main(): 
    f = makeFuncs() 
 
    print "A00(0,0):", ambgfun(f, 0, 0, 0, 0) 
    print "A01(0,0):", ambgfun(f, 0, 1, 0, 0) 
    print "A34(0,0):", ambgfun(f, 3, 4, 0, 0) 
 
if __name__ == '__main__': 
    main() 

Hermite 函数是正交的，因此所有积分都应为零。然而，它们不是，如输出所示:

A00(0,0): (1.7724538509055159, 1.4202636805184462e-08) 
A01(0,0): (8.465450562766819e-16, 8.862237123626351e-09) 
A34(0,0): (-10.1875, 26.317246925873935) 

我怎样才能使这个计算更准确？ scipy 中的 hermite 函数包含一个权重变量，该变量应用于高斯正交，如文档 (http://docs.scipy.org/doc/scipy/reference/special.html#orthogonal-polynomials) 中所给出。但是，我没有在文档中找到如何使用这些权重的提示。

希望你能帮上忙:)

谢谢，麦克斯

请您参考如下方法:

答案是您得到的结果在数值上尽可能接近于零。我不认为如果你使用 float 真的可以获得更好的结果 --- 你面临着数值积分中的一个普遍问题。

考虑一下:

import numpy as np 
from scipy import integrate, special 
f = lambda t: np.exp(-t**2) * special.eval_hermite(12, t) * special.eval_hermite(16, t) 
 
abs_ig, abs_err = integrate.quad(lambda t: abs(f(t)), -np.inf, np.inf) 
ig, err = integrate.quad(f, -np.inf, np.inf) 
 
print ig 
# -10.203125 
print abs_ig 
# 2.22488114805e+15 
print ig / abs_ig, err / abs_ig 
# -4.58591912155e-15  1.18053770382e-14 

因此，被积函数的值的计算精度与浮点 epsilon 相当。由于在减去大幅值振荡被积函数的值时存在舍入误差，因此不太可能获得更好的结果。

那么如何进行呢？根据我的经验，您现在需要做的不是从数字上而是从分析上解决问题。重要的是，厄米多项式乘以权重函数的傅立叶变换是已知的，因此您可以在这里一直在傅立叶空间中工作。