使用來自scipy.integrate的四元數值積分結果中的虛假陷波

我試圖從數字上解決使用NumPy和quad從scipy.integrate以下積分。該代碼是有點兒工作，但我得到虛假的缺口中得到的數據：使用來自scipy.integrate的四元數值積分結果中的虛假陷波

Output

任何人有他們爲什麼發生，以及如何得到正確的結果，平穩的想法？

這裏是在Python原代碼：

#!/usr/bin/env python 

import numpy as np 
from scipy.integrate import quad 
import matplotlib.pyplot as pl 

numpts = 300 

t_min = -4 
t_max = 100 
tt = np.linspace(t_min, t_max, numpts) 

mean = 0.   # ps 
fwhm = .05   # ps 

def gaussian(x, mean, fwhm): 
    return 1./np.sqrt(np.pi)/fwhm * np.exp(-1. * (x - mean)**2/fwhm**2) 

def integrand(t_, t, mean, fwhm): 
    denum = np.sqrt(t - t_) 
    r = gaussian(t_, mean, fwhm)/denum 
    return r 

def integrate(t, mean, fwhm, tmin): 
    return quad(integrand, tmin, t - 1e-9, args=(t, mean, fwhm))[0] 

if __name__ == '__main__': 
    vec_int = np.vectorize(integrate) 
    y = vec_int(tt, mean, fwhm, tt.min()) 

    fig = pl.figure() 
    ax = fig.add_subplot(111) 
    ax.plot(tt, y, 'bo-', mec='none') 
    ax.set_xlim(-5, 101) 
    pl.show()

來源

2017-02-17 François

你嘗試，並期待在額外的信息四回報？您可以將您的集成函數添加到全局（例如列表）中，然後在向量循環之後檢查它。 –

所以我通過使用tanh-sinh方法切換到mpmath庫和它自己的集成quad來解決我的問題。我也必須拆分積分，以便數據是單調的。輸出看起來是這樣的：

Output

我不是100％肯定，爲什麼這個工作，但它可能與數字精度和集成方法的行爲做。

這裏的工作代碼：

#!/usr/bin/env python 

import numpy as np 
import matplotlib.pyplot as pl 
import mpmath as mp 

numpts = 3000 

t_min = -4 
t_max = 100 
tt = np.linspace(t_min, t_max, numpts) 

mean = mp.mpf('0')   # ps 
fwhm = mp.mpf('.05')   # ps 

def gaussian(x, mean, fwhm): 
    return 1./mp.sqrt(np.pi)/fwhm * mp.exp(-1. * (x - mean)**2/fwhm**2) 

def integrand(t_, t, mean, fwhm): 
    denum = np.sqrt(t - t_) 
    g = gaussian(t_, mean, fwhm) 
    r = g/denum 
    return r 

def integrate(t, mean, fwhm, tmin): 
    t = mp.mpf(t) 
    tmin = mp.mpf(tmin) 
    # split the integral because it can only handle smooth functions 
    res = mp.quad(lambda t_: integrand(t_, t, mean, fwhm), 
        [tmin, mean], 
        method='tanh-sinh') + \ 
       mp.quad(lambda t_: integrand(t_, t, mean, fwhm), 
         [mean, t], 
         method='tanh-sinh') 
    ans = res 
    return ans 

if __name__ == '__main__': 
    mp.mp.dps = 15 
    mp.mp.pretty = True 
    vec_int = np.vectorize(integrate) 
    y = vec_int(tt, mean, fwhm, tt.min()) 

    fig = pl.figure() 
    ax = fig.add_subplot(111) 
    # make sure we plot the real part 
    ax.plot(tt, map(mp.re, y), 'bo-', mec='none') 
    ax.set_xlim(-5, 101) 
    pl.show()

來源

2017-03-07 00:39:06

我懷疑將是一個奇點積被搞亂了的四（包）的內部工作。然後我會嘗試（按此順序）：在quad中使用weights =「cauchy」;添加和減去奇點;看看告訴quadpack的方法，積分有一個難點。

來源

2017-02-17 05:57:44

使用來自scipy.integrate的四元數值積分結果中的虛假陷波

回答

相關問題