使用来自scipy.integrate的四元数值积分结果中的虚假陷波

Question

我试图从数字上解决使用NumPy和quad从scipy.integrate以下积分。该代码是有点儿工作，但我得到虚假的缺口中得到的数据：

任何人有他们为什么发生，以及如何得到正确的结果，平稳的想法？

这里是在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() 

Answer 1

所以我通过使用tanh-sinh方法切换到mpmath库和它自己的集成quad来解决我的问题。我也必须拆分积分，以便数据是单调的。输出看起来是这样的：

我不是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() 

Answer 2

我怀疑将是一个奇点积被搞乱了的四（包）的内部工作。然后我会尝试（按此顺序）：在quad中使用weights =“cauchy”;添加和减去奇点;看看告诉quadpack的方法，积分有一个难点。