我已经使用Area、sigmax和sigmay参数定义了2D高斯(自变量之间没有相关性)。当我在两个变量中从(-inf,inf)进行积分时,我只在sigmax和sigmay为1时得到Area。
import numpy as np
import scipy.integrate as sci
class BeamDistribution(object):
def __init__(self, Ipeak, sigmax, sigmay):
print Ipeak, sigmax, sigmay
self.__Ipeak = Ipeak
self.__sigmax = sigmax
self.__sigmay = sigmay
def value(self, x, y):
factor = self.__Ipeak/(2.*np.pi*self.__sigmax * self.__sigmay)
factorx = np.exp(-x**2/(2.*self.__sigmax**2))
factory = np.exp(-y**2/(2.*self.__sigmay**2))
return factor*factorx*factory
def integral(self, a, b, c, d):
integration = sci.dblquad(self.value, a, b, lambda x: c, lambda x: d,
epsrel = 1e-9, epsabs = 0)
# sci.quad_explain()
return integration
def __call__(self, x, y):
return self.value(x, y)
if __name__ == "__main__":
Ipeak = 65.0e-3
sigmax = 0.2e-3
sigmay = 0.3e-3
limit = np.inf
my_beam_class = BeamDistribution(Ipeak, sigmax, sigmay)
total = my_beam_class.integral(-limit, limit, -limit, limit)
print "Integrated total current ",total," of Ipeak ", Ipeak
my_beam_class = BeamDistribution(Ipeak, 1, 1)
total = my_beam_class.integral(-limit, limit, -limit, limit)
print "Integrated total current ",total," of Ipeak ", Ipeak
输出为
0.065 0.0002 0.0003
Integrated total current (7.452488478001055e-32, 6.855160478762106e-41) of Ipeak 0.065
0.065 1 1
Integrated total current (0.4084070449667172, 1.0138233535120856e-11) of Ipeak 0.065
知道为什么会发生这种事吗?我想这应该是一些简单的东西,但看了几个小时后,我看不出有什么错。
积分高斯核的正确方法是使用高斯-埃尔米特象限,请参见此处
它在Python中实现为SciPy模块http://docs.scipy.org/doc/numpy-1.10.1/reference/generated/numpy.polynomial.hermite.hermgauss.html
代码,高斯核上的积分应该是;(π)
import math
import numpy as np
a,w = np.polynomial.hermite.hermgauss(32)
print(a)
print(w)
def f(x):
return 1.0
s = 0.0
for k in range(0,len(a)):
s += w[k]*f(a[k])
print(s - math.sqrt(math.pi))
2D案例
Ipeak = 0.065
sigmax = 0.2e-3
sigmay = 0.3e-3
sqrt2 = math.sqrt(2.)
def h(x, y):
return Ipeak*1.0/math.pi
s = 0.0
for k in range(0, len(a)):
x = sqrt2 * sigmax * a[k]
t = 0.0
for l in range(0, len(a)):
y = sqrt2 * sigmay * a[l]
t += w[l] * h(x, y)
s += w[k]*t
print(s)
我认为0.002西格玛的高斯对于正交来说太峰值了:Scipy忽略了这个很小的峰值,到处只看到零。您有两个解决方案:
-
重新规范化函数:∫ab
(x)dx=&西格玛&inta/σb/&西格玛f(u)du -
把积分切成许多块。这里有一个例子,它计算从-无穷大到-4*西格玛,然后从-4*西格玛到4*西格玛,再从4*西格玛到无穷大的积分:
def integral(self, a, b, c, d):
integration =0
nsigmas=4
for intervalx in [(a,-nsigmas*sigmax),(-nsigmas*sigmax,nsigmas*sigmax),(nsigmas*sigmax,b)]:
for intervaly in [(c,-nsigmas*sigmay),(-nsigmas*sigmay,nsigmas*sigmay),(nsigmas*sigmay,d)]:
integration+= sci.dblquad(self.value, intervalx[0], intervalx[1], lambda x: intervaly[0], lambda x: intervaly[1],
epsrel = 1e-9, epsabs = 0)[0]
# sci.quad_explain()
return integration
我得到这个输出:
0.065 0.0002 0.0003
Integrated total current 0.06499999987174367 of Ipeak 0.065
0.065 1 1
Integrated total current 0.06500000000019715 of Ipeak 0.065