使用来自 scipy.integrate 的四边形数值积分结果中的虚假缺口

Question

我正在尝试使用 NumPy 和 scipy.integrate 中的 quad 对以下积分进行数值求解。代码有点工作，但我在结果数据中得到了虚假的缺口：

有人知道它们为什么会发生以及如何获得正确的平滑结果吗？

这里是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

我怀疑可积奇点正在扰乱 quad(pack) 的内部运作。然后我会尝试（按此顺序）：在 quad 中使用 weights="cauchy"；加减奇点；看看告诉quadpack积分有难点的方法

Answer 2

所以我通过使用 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()

使用来自 scipy.integrate 的四边形数值积分结果中的虚假缺口

Spurious notches in numerical integration results using quad from scipy.integrate

python

numpy

scipy

numerical-methods

numerical-integration