1

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:

enter image description here

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)
Dipole
  • 873
  • 2
  • 10
  • 23
  • 1
    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
  • @Kirill Please see above, yes the other question was me trying to understand Gauss-Laguerre and how I could apply it to this problem. I'm not sure how to go about researching numerical techniques because all I know is that the intregand in question can be described as "singular" and "oscillatory." – Dipole Sep 23 '14 at 21:23
  • @Jack Right, I missed that. But in general I expect the integration routines themselves to handle singularities, not to do it myself. Also, another source of error might be the linear interpolation in the integrand. – Kirill Sep 23 '14 at 21:34
  • @Kirill Yes I was thinking that too. I made sure that the integrand was sampled very finely when I calculated it so I assumed linear interpolation would have been fine, but perhaps I should use a higher order interpolation, maybe cubic? – Dipole Sep 23 '14 at 21:36
  • The problem with linear/cubic interpolation might be that any linear/cubic interpolant has no singularities at all, so it wouldn't approximate the integrand in a fundamental way. Is it at all possible to sample the integrand at will? (Also, stackexchange is not really such a great platform for long-winded questions with lots of discussion; it's really mainly designed for questions with unique clear answers.) – Kirill Sep 23 '14 at 21:39
  • @Kirill I think I see what you mean, but the sampled values I have for functions involved in the integrand are so finely sampled, and that is what is interpolated. The singularity comes from the $\frac{1}{\sqrt{\zeta^2(s)-1}}$ which I of course call directly each time, only with $\zeta(s)$ being interpolated. Hence the singularity should be preserved right? – Dipole Sep 23 '14 at 21:49
  • Yes, if the interpolant achieves the value $1$ at $s=0$. – Kirill Sep 23 '14 at 21:52
  • @Kirill Yes the value of $\zeta = 1$ at $s=0$ because the trajectory of $\zeta(s)$ is calculated by starting from there and stepping forward in $s$, then if I need intermediate points I just interpolate. – Dipole Sep 23 '14 at 21:56
  • The relative error of $\zeta(s)^2-1$ when you know $\zeta(s)^2$ to within $\pm\delta$ is $\delta (\zeta(s)^2-1)^{-1}$. So the absolute error of $(\zeta(s)^2-1)^{-1/2}$ will be $(\zeta(s)^2-1)^{-3/2}\delta/2$, which is large, singular, and not integrable. I think this is the main issue with interpolation here, the integrand is (more) ill-conditioned. – Kirill Sep 23 '14 at 22:05
  • @Kirill Ok so I guess I will have to modify my code so that I call the function which calculates $\zeta(s)$, $t_1(s)$, $t_2(s)$ on every call by quad. Is there any other fault you can see otherwise with my implementation above? Many thanks, I really appreciate your help! – Dipole Sep 23 '14 at 22:12
  • Try to get rid of the cancellation in $\zeta(s)^2-1$ by finding some other way to compute it that doesn't involve subtraction of nearly-equal approximately-known quantities. Use something standard, like github gists for code/data downloads. Let the integration routine handle the singularities itself, somebody already wrote it and gave it some thought. There are too many potential issues here for someone to figure out exactly what's going on, let alone give a clean definitive answer. – Kirill Sep 23 '14 at 22:24
  • @Kirill But is it even possible to get rid of this cancellation? I take it you mean something of the like discussed in: http://math.stackexchange.com/questions/136634/rewriting-to-avoid-catastrophic-cancellation. I guess I could start with something similar and factorize my expression. – Dipole Sep 23 '14 at 22:34
  • @Kirill ok, ill brainstorm, thanks. But just to clarify, I can't rewrite the definition of $\zeta$ as it is computed numerically. – Dipole Sep 23 '14 at 22:46
  • If you can't rewrite ζ to get rid of cancellation, it's just an ill-conditioned problem, so probably not much you can do. OTOH, "maybe" the integral doesn't depend that much on that value of $\zeta$ and it might still be accurate. – Kirill Sep 23 '14 at 22:52
  • The anwer in your linked question is, I think, misleading because it picks some convenient numbers. If $x, y$ are approximate, the issue is that both $x^2-y^2$ and $x-y$ are inaccurate; that answer assumes that both $x,y$ are also very large, which magnifies the inaccuracy of the first expression. In your case that's not relevant, and writing $\zeta^2-1=(\zeta-1)(\zeta+1)$ does nothing at all. – Kirill Sep 23 '14 at 22:56
  • @Kirill ok, I guess that is the case then, because $\zeta$ was calculated by numerical integration of an ODE which has no analytical solution. Unfortunately scipy.quad does not like when I set the integration interval to $(0,\infty)$... – Dipole Sep 23 '14 at 22:58
  • 2
    @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 01:29
  • Comments are not for extended discussion; this conversation has been moved to chat. – Paul Sep 24 '14 at 01:55
  • @Kirill I Updated my question. I was able to find an analytical expression for my integrand. – Dipole Oct 09 '14 at 10:45

0 Answers0