Falsche Kerben in numerischen Integrationsergebnissen unter Verwendung von Quad von scipy.integrate

Question

Ich versuche numerisch das folgende Integral mit NumPy und quad von scipy.integrate zu lösen. Der Code ist ein bisschen arbeiten, aber ich unechte Kerben in den resultierenden Daten:

Jeder hat eine Idee, warum sie auftreten und wie man das richtige glatte Ergebnis zu bekommen?

Hier ist der ursprüngliche Code in 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

Also löste ich mein Problem, indem ich auf die mpmath Bibliothek und seine eigene Integration quad mit der tanh-sinh Methode umschaltete. Ich musste auch das Integral aufteilen, so dass die Daten monoton sind. Die Ausgabe sieht wie folgt aus:

Ich bin nicht 100% sicher, warum dies funktioniert, aber es kann mit der numerischen Genauigkeit und das Verhalten des Integrationsmethode zu tun haben.

Hier ist der Arbeitscode:

#!/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

Mein Verdacht wäre, dass eine integrierbare Singularität ist das Innenleben von Quad (Pack) vermasselt. Ich würde dann versuchen (in dieser Reihenfolge): Verwenden Sie Gewichte = "Cauchy" in Quad; Addiere und subtrahiere die Singularität; schau dir an, wie man Quadpack sagen kann, dass das Integral einen schwierigen Punkt hat.