神经质的集成.quad未返回预期值

Question

我有下面的函数，我想用python进行数值积分

使用scipy，我编写了以下代码：

def voigt(a,u):

fi = 1
er = Cerfc(a)*np.exp(np.square(a))
c1 = np.exp(-np.square(u))*np.cos(2*a*u)
c1 = c1*er  #first constant term
pis = np.sqrt(np.pi) 

c2 = 2./pis    #second constant term    

integ  = inter.quad(lambda x: np.exp(-(np.square(u)-   
 np.square(x)))*np.sin(2*a*(u-x)), 0, u)

print integ
ing = c1+c2*integ[0]

return ing

对于Cerfc（a）函数，我只使用scipy.erfc公司计算互补误差函数。

所以这个函数对于u的值很小，但是大的u值（超过60 ish）会破坏代码，并且会得到非常小的数字。例如，如果输入a=0.01和u=200，结果是1.134335928072937e-40，其中真实答案是：1.410526851411200e−007

除此之外，四元计算的错误scipy返回与答案的顺序相似。我真的很为难你，我真的很感激你的帮助。

这是家庭作业，但这是物理课。所以这个计算只是物理学中更广泛问题的一个步骤。如果你帮我，你就不会帮我作弊了：）

Answer 1

根据维基百科文章Voigt profile，Voigt functions U（x，t）和V（x，t）可以用络合物Faddeeva functionw（z）表示：

U(x,t) + i*V(x,t) = sqrt(pi/(4*t))*w(i*z)

Voigt函数H（a，u）可以用u（x，t）表示为

^{pr2}$

（另请参见DLMF section on Voigt functions。）

scipy在^{}中有一个Faddeeva函数的实现。 下面是Voigt函数的实现：

from __future__ import division

import numpy as np
from scipy.special import wofz


_SQRTPI = np.sqrt(np.pi)
_SQRTPI2 = _SQRTPI/2

def voigtuv(x, t):
    """
    Voigt functions U(x,t) and V(x,t).

    The return value is U(x,t) + 1j*V(x,t).
    """
    sqrtt = np.sqrt(t)
    z = (1j + x)/(2*sqrtt)                    
    w = wofz(z) * _SQRTPI2 / sqrtt
    return w

def voigth(a, u):
    """
    Voigt function H(a, u).
    """
    x = u/a
    t = 1/(4*a**2)
    voigtU = voigtuv(x, t).real
    h = voigtU/(a*_SQRTPI)
    return h

你说过当a=0.01，u=200时，H（a，u）的值是1.410526851411200e−007。我们可以检查：

In [109]: voigth(0.01, 200)
Out[109]: 1.41052685142231e-07

---

上面没有回答当u很大时，代码为什么不能工作的问题。要成功地使用quad，最好对被积函数有一个很好的理解。在您的例子中，当u很大时，只有x = u附近的一个很小的区间对积分有很大的贡献。quad没有检测到这一点，因此它忽略了整数的大部分，并返回一个太小的值。在

解决此问题的一种方法是使用quad的points参数，该点非常接近间隔的终点。例如，我将quad的调用改为：

integ = inter.quad(lambda x: np.exp(-(np.square(u)-np.square(x))) * np.sin(2*a*(u-x)),
                   0, u, points=[0.999*u])

在这种变化下，函数返回voigt(0.01, 200)：

In [191]: voigt(0.01, 200)
Out[191]: 1.4105268514252487e-07

对于值0.999*u，我没有一个严格的理由；这仅仅是一个距离间隔末尾足够近的点，可以对大约200左右的u给出一个合理的答案。进一步研究被积函数可以给你一个更好的选择。（例如，你能找到被积函数最大值位置的解析表达式吗？如果是这样，那将比0.999*u好得多

您也可以尝试调整epsabs和epsrel的值，但是在我的几个实验中，添加points参数产生了最大的影响。在