scipy "The integral is probably divergent, or slowly convergent"

Question

R = 0.05; l = 0.1;
q0 = 200; 
alpha = 17.64;
lambda1 = 67.9;

def integral2(x, z, r):
  s = 0;
  if z > 0:
    s = 1
  elif z == 0:
    s = 0.5
  else: 
    s = 0
  B = ((lambda1*math.tan(x) - alpha) * math.exp(2 * math.tan(x) * l)) - ((lambda1*math.tan(x) + alpha) * math.exp((-2) * math.tan(x) * l))
  BZ = ((math.cosh(z+l) * ((lambda1 * math.tan(x) * math.cosh(math.tan(x)*l)) - (alpha * math.sinh(math.tan(x)*l)))) / B) - ((math.sinh(math.tan(x)*z) * s) / math.tan(x))
  result = ((mpmath.j0(r * math.tan(x)) * mpmath.j1(R * math.tan(x))) / (math.pow(math.cos(x),2) * math.pow(math.tan(x),2))) * BZ
  return result

r = 0
resArr = [];
rArr = [];
while r < 1:
  r += 0.01
  r = math.ceil(r*100)/100;
  res, _ = scipy.integrate.quad(integral2, 0, (np.pi / 2) - 0.02741184711388415, args=(l/2, r))
  resArr.append(res*R*q0/lambda1)
  rArr.append(r)

Hello! I am trying to understand why I get "The integral is probably divergent, or slowly convergent." in this case. Can anybody help me?

Thanks!

Answer 1

After adding the required imports, we can plot your integrand:

import numpy as np
import mpmath
import math
import matplotlib.pyplot as plt
import scipy.integrate

R = 0.05; l = 0.1;
q0 = 200; 
alpha = 17.64;
lambda1 = 67.9;

def integral2(x, z, r):
  s = 0;
  if z > 0:
    s = 1
  elif z == 0:
    s = 0.5
  else: 
    s = 0
  B = ((lambda1*math.tan(x) - alpha) * math.exp(2 * math.tan(x) * l)) - ((lambda1*math.tan(x) + alpha) * math.exp((-2) * math.tan(x) * l))
  BZ = ((math.cosh(z+l) * ((lambda1 * math.tan(x) * math.cosh(math.tan(x)*l)) - (alpha * math.sinh(math.tan(x)*l)))) / B) - ((math.sinh(math.tan(x)*z) * s) / math.tan(x))
  result = ((mpmath.j0(r * math.tan(x)) * mpmath.j1(R * math.tan(x))) / (math.pow(math.cos(x),2) * math.pow(math.tan(x),2))) * BZ
  return result

r = 0.01
x = np.linspace(1e-6, (np.pi / 2) - 0.02741184711388415, 300)
y = [integral2(xi, l/2, r) for xi in x]
plt.plot(x, y)

At x=0, there is a singularity.

If we zoom in a bit, we can see there appears to be another singularity around x=0.85.

Even if they are integrable, they are going to be challenging for numerical integrators. This is why you're getting the warning: quad knows its having trouble, and it's letting you know that the results may not be reliable. For scipy.integrate.quad to have any hope of handling the mid-interval singularity, you'd need to use the points argument to tell it where the singularity occurs. I decided to go a different route to investigate further, though.

Some quadrature schemes can handle some singularities at the endpoints of the integration interval. mpmath's implementation of the Tanh-Sinh rule is one such scheme, and since you were already using mpmath (albeit unnecessarily - scipy.special has j0 and j1), I figured I'd see if mpmath.quad could handle the singularity at x=0.

from mpmath import mp
mp.dps = 50

R = mp.mpf(0.05)
l = mp.mpf(0.1)
q0 = mp.mpf(200.)
alpha = mp.mpf(17.64)
lambda1 = mp.mpf(67.9)

def integrand(x, z, r):
  s = mp.mpf(0.);
  if z > 0:
    s = mp.mpf(1.)
  elif z == 0:
    s = mp.mpf(0.5)

  B = ((lambda1*mp.tan(x) - alpha) * mp.exp(2 * mp.tan(x) * l)) - ((lambda1*mp.tan(x) + alpha) * mp.exp((-2) * mp.tan(x) * l))
  BZ = ((mp.cosh(z+l) * ((lambda1 * mp.tan(x) * mp.cosh(mp.tan(x)*l)) - (alpha * mp.sinh(mp.tan(x)*l)))) / B) - ((mp.sinh(mp.tan(x)*z) * s) / mp.tan(x))
  result = ((mp.j0(r * mp.tan(x)) * mp.j1(R * mp.tan(x))) / (mp.cos(x)**2 * mp.tan(x)**2)) * BZ
  return result

r = mp.mpf(0.01)
l = mp.mpf(l)
# integrate just from 0 to 0.8
res = mp.quad(lambda x: integrand(x, l/2, r), (0, 0.8), method='tanh-sinh')
print(res)

The result appears to double when the working precision is doubled, so I don't think it is converging to a correct answer. Your integral might not be finite.