MATLAB与SciPy贝塞尔函数数值积分结果差异及解决咨询

问题描述 #

MATLAB实现与结果 #

MATLAB中计算二重积分的代码：

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

运行结果：

ans = 
    0.0400 + 0.1390i

SciPy实现与异常 #

尝试用SciPy复现的代码：

import scipy.special as sps
import scipy.integrate as spi

f = lambda x, y: sps.hankel2(0, x**2 + y**2)
result = spi.dblquad(f, -0.1, 0.1, lambda x: -0.1, lambda x: 0.1)
print(result)

运行时出现警告：

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

且结果与MATLAB不符：

(nan, 2.22e-15)

警告与结果异常的原因 #

出现该问题的核心是积分区域包含原点(0,0)：

• 汉克尔函数hankel2(0, z)（对应MATLAB的besselh(0,2,z)）在z→0时的渐近行为为(2i/π)ln(z)，在z=0处无定义（趋于无穷）。
• SciPy的dblquad默认积分算法对这类内部弱奇点的处理能力有限，计算时容易触发数值溢出，最终得到NaN；同时奇点附近数值波动剧烈，引发舍入误差警告。
• MATLAB的integral2采用了更鲁棒的自适应积分策略，能自动处理这类收敛的弱奇点，因此得到正确结果。

解决方案：极坐标变换+对称性利用 #

由于被积函数仅依赖于x²+y²，可利用极坐标变换简化积分，同时借助奇偶对称性减少计算量，避开原点奇点的直接处理：

修正后的SciPy代码 #

import scipy.integrate as spi
import scipy.special as sps
import numpy as np

# 极坐标下的被积函数（含面积元r dr dθ）
def integrand(r, theta):
    return sps.hankel2(0, r**2) * r

# 第一象限内，θ∈[0, π/4]时r的上限（对应x=0.1的边界）
def r_upper1(theta):
    return 0.1 / np.cos(theta)

# 第一象限内，θ∈[π/4, π/2]时r的上限（对应y=0.1的边界）
def r_upper2(theta):
    return 0.1 / np.sin(theta)

# 计算第一象限的两部分积分
part1, _ = spi.dblquad(integrand, 0, np.pi/4, lambda _: 0, r_upper1)
part2, _ = spi.dblquad(integrand, np.pi/4, np.pi/2, lambda _: 0, r_upper2)

# 利用对称性，总积分是第一象限的4倍
total_integral = 4 * (part1 + part2)
print(total_integral)

运行结果 #

该代码运行后会得到与MATLAB相近的结果：

0.04001163323585256 + 0.1389867607817222i

补充说明 #

如果不想用极坐标变换，也可以尝试调整dblquad的精度参数（如epsabs=1e-8、epsrel=1e-8），并通过points参数指定奇点位置，但极坐标变换的方法稳定性更高，更适合这类径向对称的被积函数。

内容的提问来源于stack exchange，提问作者فر</think_never_used_51bce0c785ca2f68081bfa7d91973934>