Локальные вырезы в результатах цифровой интеграции с использованием quad от 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

Итак, я решил свою проблему, переключившись на библиотеку mpmath и свою собственную интеграцию quad с использованием метода tanh-sinh. Мне также пришлось разбить интеграл так, чтобы данные были монотонны. Результат выглядит следующим образом:

Я не 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

Подозреваю бы, что интегрируемая особенность портя внутреннюю работу четырехугольника (пакет). Затем я попытался (в этом порядке): использовать weight = «cauchy» в quad; добавить и вычесть сингулярность; посмотрите на способы сказать квадпаку, что интеграл имеет сложную точку.