# Using scipy.quad to calculate difficult integral

When evaluating the integral below in python using scipy.quad I get the following warning:

UserWarning: The maximum number of subdivisions (50) has been achieved. If increasing the limit yields no improvement it is advised to analyze the integrand in order to determine the difficulties. If the position of a local difficulty can be determined (singularity, discontinuity) one will probably gain from splitting up the interval and calling the integrator on the subranges. Perhaps a special-purpose integrator should be used.

Furthermore the real integration bounds should be zero to infinity (see below) but when I change the bounds to this my result looks very different from the finite values used in my example below. How should I go about calculating this integral?

I have the following complex integral to compute numerically:

\begin{align} \int_0^{\infty} Re \left( \frac{t_1(s)-it_2(s)}{\sqrt{\zeta^2(s)-1}} \right) \exp(i\omega s/V) \, \mathrm{d}s \end{align}

here $i=\sqrt{-1}$, $t_1$ and $t_2$ are real parametric functions of $s$, $V=\frac{\pi}{2(\pi+2)}$ is a constant, and $\zeta(s)$ is a complex function of $s$. The functions $t_1$, $t_2$, and $\zeta$ are not available in closed form, rather I have computed them numerically. A numpy array containing the sampled function values at various vales of s can be found to download at http://speedy.sh/prG9e/stackexchange.npy. The integrand is problematic because $\zeta(0) = 1$.

I have plotted the integrands for various values of omega:

The singularity in the real integrand is apparent. My work so far:

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate

T2 = np.load('stackexchange.npy') #Data: [s, x1, x2, t1, t2]

def pyfunc(z):
return np.sqrt(z**2-1)

def integrand(s, omega):
x1 = np.interp(s, T2[:,0], T2[:,1] )
x2 = np.interp(s, T2[:,0], T2[:,2] )
t1 = np.interp(s, T2[:,0], T2[:,3] )
t2 = np.interp(s, T2[:,0], T2[:,4] )
zeta = x2+x1*1j
sigma = np.pi/(np.pi+2)
V = 1/(2*sigma)
return (-t2*np.real(1j/pyfunc(zeta))+t1*np.real(1/pyfunc(zeta)))*np.exp(1j*omega*s/V)

def integral(omega):
def real_func(x,omega):
return np.real(integrand(x,omega))
def imag_func(x,omega):
return np.imag(integrand(x,omega))
a = 0.05 #Lower bound
b = 20.0 #Upper bound
real_integral = integrate.quad(real_func, a, b, args=(omega))
imag_integral = integrate.quad(imag_func, a, b, args=(omega))
return real_integral[0] + 1j*imag_integral[0]

vintegral = np.vectorize(integral)

omega = np.linspace(-30, 30, 601)
I = integral(omega)

plt.plot(omega, I.real, omega, I.imag)
plt.show()


# Edit

I found an analytical representation of the integrand in question, defined below. It still seems that there are some numerical difficulties with it though as NaNs are returned when I try to integrate.

    def pyfunc(z):
return np.sqrt(z**2-1)

def func(theta):
t1 = 1/np.sqrt(1+np.tan(theta)**2)
t2 = -1/np.sqrt(1+1/np.tan(theta)**2)
return t1, t2

def integrand(s, omega):
sigma = np.pi/(np.pi+2)
xs = np.exp(-np.pi*s/(2*sigma))
x1 = -2*sigma/np.pi*(np.log(xs/(1+np.sqrt(1-xs**2)))+ np.sqrt(1-xs**2))
x2 = 1-2*sigma/np.pi*(1-xs)
zeta = x2+x1*1j
Vc = 1/(2*sigma)
theta =  -1*np.arcsin(np.exp(-np.pi/(2.0*sigma)*s))
t1, t2 = func(theta)
return np.real((t1-1j*t2)/pyfunc(zeta))*np.exp(1j*omega*s/Vc)

def integral(omega):
def real_func(x,omega):
return np.real(integrand(x,omega))
def imag_func(x,omega):
return np.imag(integrand(x,omega))
a = 0
b = np.inf
real_integral = integrate.quad(real_func, a, b, args=(omega))
imag_integral = integrate.quad(imag_func, a, b, args=(omega))
return real_integral[0] + 1j*imag_integral[0]

vintegral = np.vectorize(integral)

• You may want to look at the suggestions in Method for numerical integration of difficult oscillatory integral. – Geoff Oxberry Sep 23 '14 at 18:00
• I don't understand: why are the bounds $(0.05,20)$, not $(0,\infty)$? One thing to try, because of the singularity present, is applying trapezoid rule after the double-exponential substitution $s=\exp(c\sinh u)$. Also, depending on how close $\zeta(s)$ is to $1$ when $s$ is small, you might have big relative errors in the term $\zeta(s)^2-1$. Can you post a "similar" closed-form integrand that has the same major properties your integrand has? – Kirill Sep 23 '14 at 21:02
• By the way, in your previous question you asked about Gauss-Laguerre integration for $e^{(1-s^2)t}\,dt$; here the integrand is even worse because the $2n$-th derivative of $\frac1{\sqrt{t}}$ is super-exponential in $n$. – Kirill Sep 23 '14 at 21:09
• Perhaps I should have been a bit more clear, but I stated that when I used $(0,\infty)$ I get some problems, so I decided to first of all exclude the singularity (0.05 was just a small arbitrary number, and 20 was supposed to be high enough such that it incurs a small error. Besides I only have samples of my functions up to $s=30$ if you look at the data array. Im not sure how I could find a similar closed form integral, nor how to use the substation you have provided. I appreciate the advice though, ill continue to read up on what you suggest. – Dipole Sep 23 '14 at 21:11
• @Jack: Agreed with Kirill here: Stack Exchange isn't good for long questions with lots of discussion. Having lots of comments in the comment section will get flagged (like this one). In addition, it's generally frowned upon to post repeat questions that are variations of each other. Your posts aren't exactly the same, but they're not too far off, and posting too many more about integration in a short span will start triggering mods to close them. As for the discussion, it's better to have a discussion in a chat room than in the comments. – Geoff Oxberry Sep 24 '14 at 1:29