我想用scipy.integrate.quad. 在
然而,由于包含了大阶多项式,因此存在数值误差。我的代码是:
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函数是正交的,因此所有的积分都应该等于零。但是,它们不是,如输出所示:
^{pr2}$我怎样才能使这个计算更准确?scipy的hermite函数包含一个用于高斯求积的权值变量,如文档(http://docs.scipy.org/doc/scipy/reference/special.html#orthogonal-polynomials)中所示。然而,我在文档中没有发现如何使用这些权重的提示。在
希望您能帮忙:)
谢谢,麦克斯
答案是,你得到的结果在数值上是接近于零的。我不认为如果你用浮点数来处理,你真的不可能得到更好的结果——你在数值积分中面临着一个普遍的问题。在
考虑一下这个:
因此,被积函数的计算精度与浮点epsilon相当。由于在大幅度振荡被积函数的减法中存在舍入误差,因此不可能得到更好的结果。在
那么如何继续呢?以我的经验,你现在需要做的是解决这个问题,而不是从数字上,而是从分析上。重要的是,Hermite多项式的Fourier变换乘以权函数是已知的,所以你可以一直在Fourier空间中工作。在
相关问题 更多 >
编程相关推荐