SciPy的integrate.quad不返回

我有以下的功能，我想數字使用python集成，SciPy的integrate.quad不返回

使用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

除此之外，quad計算的錯誤scipy返回值是按照與答案類似的順序。我真的很難過，非常感謝你的幫助。

這是一個家庭作業，但它是一個物理課程。所以這個計算只是物理學中一個更廣泛問題的一個步驟。你會不會幫我作弊，如果你幫我:)

來源

2016-11-10 handroski

使用'quad'或其他特定的scipy函數是它的一部分嗎？或者您是否需要Voigt功能的實施？ –

它只是工作。謝謝 – handroski

根據維基百科的文章Voigt profile，在複雜的Faddeeva function瓦特方面Voigt functions U（X，T）和V（X，T）可以表示（Z）：

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

的Voigt函數H（A，U）可以在U（X，t）的條款

H(a,u) = U(u/a, 1/(4*a**2))/(a*sqrt(pi))

（另請參見DLMF section on Voigt functions。）

表示0

scipy在scipy.special.wofz中實現了Faddeeva函數。利用這一點，這裏的沃伊特功能的實現：

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

你說，你知道的即H值（A，U）是1.410526851411200e-007當a = 0.01和u = 200。我們可以檢查：

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值，但在我的幾個實驗中，加入了最大影響的論點。

來源

2016-11-11 00:59:41

SciPy的integrate.quad不返回

回答

相關問題