Differences Between MATLAB and SciPy numerical Integration Results with Bessel Functions

Question

I'm looking to replicate the results of the following MATLAB code in SciPy.

MATLAB Version:

f = @(x, y) besselh(0, 2, x.^2 + y.^2);
integral2(f, -0.1, 0.1, -0.1, 0.1)

The MATLAB result is:

ans = 
    0.0400 + 0.1390i

SciPy Version:

f = lambda x, y: scipy.special.hankel2(0, x**2 + y**2)
scipy.integrate.dblquad(f, -0.1, 0.1, lambda x: -0.1, lambda x: 0.1)

However, when I run the SciPy code, I receive the following warning:

IntegrationWarning: The occurrence of roundoff error is detected ...

And the result is different from the MATLAB output:

(nan, 2.22e-15)

Can someone explain why I'm getting this warning and how to achieve a result similar to MATLAB's?

Answer 1

I believe the issue is that the integral goes through the origin, which produces NaN for H0. This is because H0 = J0 - i*Y0, and Y0 asymptotes the y-axis.

For this specific case, you can use the H0 definition above to split the integrals into the real and complex parts. When doing this, you can make sure you don't cross 0 for the complex part of the integral.

import scipy

j0 = lambda x, y: scipy.special.jv(0, x**2 + y**2)
y0 = lambda x, y: scipy.special.yv(0, x**2 + y**2)
r1 = scipy.integrate.dblquad(j0, -0.1, 0.1, -0.1, 0.1)
r21 = scipy.integrate.dblquad(y0, 1e-10, 0.1, -0.1, 0.1)
r22 = scipy.integrate.dblquad(y0, -0.1, -1e-10, -0.1, 0.1)
res = r1[0] - 1j*(r21[0]+r22[0])  # 0.039999377783047595+0.13896313686573833j