Scipy：Hermite函數與正交權重的積分

Question

我想使用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() 

埃爾米特函數是正交的，因此所有的積分應該等於零。但是，他們都沒有，作爲輸出顯示：

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函數包含一個權重變量，應該用於Gaussian Quadrature，如文檔（http://docs.scipy.org/doc/scipy/reference/special.html#orthogonal-polynomials）中給出的。但是，我沒有在文檔中找到如何使用這些權重的提示。

我希望你能幫助:)

謝謝，最大

Answer 1

答案是，你得到的結果是數值接近於零，因爲它得到。我不認爲如果你使用浮點數來得到更好的結果是真的可能的 - 你在數值積分中面臨一個普遍的問題。

考慮一下：

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 

被積函數的值因此被計算爲的精度相媲美的浮點小量。由於減去大幅度振盪積分的值的舍入誤差，實際上不可能獲得更好的結果。

那麼如何繼續？根據我的經驗，你現在需要做的就是解決問題，而不是數字，但分析。重要的是，Hermite多項式乘以權函數的傅里葉變換是已知的，所以您可以在傅里葉空間中一直工作。