Scipy integration.quad未返回期望值 [英] Scipy integrate.quad not returning the expected values
问题描述
我有以下函数，我想使用python进行数字积分，
I have the following function which I would like to numerically integrate using python,
我使用scipy编写了以下代码:
Using scipy, I have written this code:
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计算互补误差函数.
For the Cerfc(a) function, I just use scipy.erfc to calculate the complimentary error function.
因此，此函数对于u的较低值确实非常有效，但是u值较大(超过60 ish)会破坏代码，并且会产生非常小的数字.例如，如果我输入a = 0.01且u = 200，则结果为1.134335928072937e-40，其中真正的答案是:1.410526851411200e−007
So this function works really well for low values of u, however larges u values (beyond 60 ish) breaks the code and I ger very very small numbers. For example, if I enter a = 0.01 and u = 200, the result is 1.134335928072937e-40, where the true answer is: 1.410526851411200e−007
除此之外，用于四边形计算的错误密码返回与答案的顺序相似.我真的很受困，非常感谢您的帮助.
In addition to this, the error scipy return for the quad calculation is on a similar order to the answer. I'm really stumped here and would really appreciate help.
这是一项家庭作业，但这是一门物理课程.因此，这一计算只是物理学中一个更广泛问题中的一步.如果您帮助我，您将不会帮助我作弊:)
This is for a homework assignment but it's a physics course. So this calculation is just one step in a broader question in physics. You will not be helping me cheat if you help me :)
推荐答案
根据Wikipedia文章 Voigt资料， Voigt函数 U(x，t)和V(x，t)可以用复数 Faddeeva函数 w(z):
According to the wikipedia article Voigt profile, the Voigt functions U(x,t) and V(x,t) may be expressed in terms of the complex Faddeeva function w(z):
U(x,t) + i*V(x,t) = sqrt(pi/(4*t))*w(i*z)
Voigt函数H(a，u)可以用U(x，t)表示为
The Voigt function H(a,u) can be expressed in terms of U(x,t) as
H(a,u) = U(u/a, 1/(4*a**2))/(a*sqrt(pi))
(另请参阅Voigt函数的 DLMF部分.)
(Also see the DLMF section on Voigt functions.)
scipy has an implementation of the Faddeeva function in scipy.special.wofz.
Using that, here's an implementation of the Voigt functions:
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.我们可以检查:
You said that you know that value of H(a,u) is 1.410526851411200e−007 when a=0.01 and u=200. We can check:
In [109]: voigth(0.01, 200)
Out[109]: 1.41052685142231e-07
上面的内容没有回答为什么u很大时为什么您的代码不起作用的问题.要成功使用quad，对被积物有很好的了解永远是一个好主意.在您的情况下，当u大时，在x = u附近只有很小的时间间隔会对积分产生重大影响. quad不会检测到它，因此它会丢失很大一部分积分并返回一个太小的值.
The above doesn't answer the question of why your code doesn't work when u is large. To use quad successfully, it is always a good idea to have a good understanding of your integrand. In your case, when u is large, only a very small interval near x = u makes a significant contribution to the integral. quad doesn't detect this, so it misses a big part of the integral and returns a value that is too small.
解决此问题的一种方法是使用quad的points自变量，其点非常接近间隔的终点.例如，我将quad的调用更改为:
One way to fix this is to use the points argument of quad with a point that is very close to the end point of the interval. For example, I changed the call of quad to:
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)返回的内容:
With that change, here's what your function returns for voigt(0.01, 200):
In [191]: voigt(0.01, 200)
Out[191]: 1.4105268514252487e-07
对于值0.999*u，我没有严格的理由；这仅是距离区间末尾足够近的一个点，可以为200左右的u提供合理的答案.对被积物的进一步调查可以为您提供更好的选择. (例如，您可以找到一个关于被积数最大值位置的解析表达式吗?如果是这样，那将比0.999*u好得多.)
I don't have a rigorous justification for the value 0.999*u; that is just a point close enough to the end of the interval to give a reasonable answer for u around 200 or so. Further investigation of the integrand could give you a better choice. (For example, can you find an analytical expression for the location of the maximum of the integrand? If so, that would be much better than 0.999*u.)
您也可以尝试调整epsabs和epsrel的值，但是在我的一些实验中，添加points参数产生了最大的影响.
You could also try tweaking the values of epsabs and epsrel, but in my few experiments, adding the points argument made the biggest impact.
这篇关于Scipy integration.quad未返回期望值的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
