Scipy積分出錯

Question

我試圖集成一個函數。該功能是保證非負：

def function(x): 
    something = ... 
    something_else = ... 
    return exp(something)/sqrt(something_else) 

現在我正在整合它：

def integrand(a, b): 
    return scipy.integrate.quad(function, a, b) 

我得到的結果不是我所期望的，所以要檢查我這樣做：

for x in range(0, 10000): 
    if integrand(0,x+1) < integrand(0,x): 
     raise ValueError("Weird!") 

果然，我越來越'奇怪'的例外。怎麼可能？

Answer 1

沒有一個工作的例子，很難調試你的代碼，但這裏有一些基於你所顯示的建議。

quad返回一個元組，其中包含兩個值：積分估計值和結果絕對誤差估計值。您似乎沒有使用代碼中的第一個值。嘗試調用quad後更改integrand到

def integrand(a, b): 
    return scipy.integrate.quad(function, a, b)[0] 

注意[0]。

接下來，當您發現「奇怪」的情況時，您應該打印integrand(0, x)和integrand(0, x+1)的值。如果它們大致相同，那麼正常的數值不精確可能會破壞您理論上預期的順序。看看quad返回的第二個值。如果該值大於integrand(0, x)和integrand(0, x+1)之間的差值，則不能真正期望積分的估計值單調遞增，因爲預期的增加量小於計算中的誤差。如果這是問題，您可以嘗試將epsabs和/或epsrel設置爲比1.49e-8的默認值更小的值，但這隻會讓您感到滿意。最終你會遇到浮點計算的限制。