使用 scipy.quad 在 python 中评估以下积分时，我收到以下警告：

用户警告：已达到最大细分数 (50)。如果增加限制没有改善，建议分析被积函数以确定困难。如果可以确定局部难度的位置（奇点、不连续性），则可能会从拆分区间并调用子范围上的积分器中获益。也许应该使用专用的积分器。

此外，真正的积分界限应该是零到无穷大（见下文），但是当我将界限更改为此时，我的结果看起来与下面示例中使用的有限值非常不同。我应该如何计算这个积分？

我有以下复杂的积分来进行数值计算：

$$\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}$$

这里$i=\sqrt{-1}$,$t_1$和$t_2$是实参数函数$s$,$V=\frac{\pi}{2(\pi+2)}$是一个常数，并且$\zeta(s)$是一个复杂的函数$s$. 功能$t_1$,$t_2$， 和$\zeta$不是以封闭形式提供的，而是我用数字计算了它们。可以在http://speedy.sh/prG9e/stackexchange.npy下载一个 numpy 数组，其中包含在 s 的不同值处采样的函数值。被积函数是有问题的，因为$\zeta(0) = 1$.

我已经绘制了各种欧米茄值的被积函数：

实被积函数的奇点是显而易见的。到目前为止我的工作：

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()

编辑

我找到了所讨论的被积函数的分析表示，定义如下。尽管在我尝试积分时返回了 NaN，但它似乎仍然存在一些数值困难。

    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)