CODEWITHSUNDEEP
pythonnumpyscipynumerical-methodsnumerical-integration
François
François

Reputation: 133

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

I am trying to solve numerically the following integral using NumPy and quad from scipy.integrate. The code is kinda working, but I get spurious notches in the resulting data:

Output

Anybody has an idea why are they occurring and how to get the correct smooth result?

Here is the original 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()

Upvotes: 0

Views: 177

Answers (2)

François
François

Reputation: 133

So I solved my problem by switching to the mpmath library and its own integration quad using the tanh-sinh method. I also had to split the integral so that the data is monotonous. The output looks like this:

Output

I am not 100% sure why this works, but it may have to do with the numeric precision and the behaviour of the integration method.

Here's the working code:

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

Upvotes: 0

ev-br
ev-br

Reputation: 26090

My suspicion would be that an integrable singularity is messing up the inner workings of quad(pack). I'd then try (in this order): use weights="cauchy" in quad; add and subtract the singularity; look at ways of telling quadpack that the integral has a difficult point.

Upvotes: 1

Related Questions