Scipy.integrate中Romberg方法在积分限为-pi/2时的异常问题

问题重现 #

使用Anaconda环境下的Scipy 1.10.0运行以下代码：

from math import cos, pi
from scipy.integrate import romberg

f = lambda x: x**2*cos(x)**2

res = romberg(f, -pi/2, pi/2)
print(res)

res = romberg(f, 0, pi/2)
print(res)

dx = 1e-4
res = romberg(f, -pi/2+dx, pi/2)
print(res)

得到输出：

9.687909744833307e-33
0.25326501581059374
0.5065300316142199

而预期的正确结果（积分∫_{-π/2}^{π/2} x²cos²x dx）应为0.5065300316211875。可见当积分下限为-pi/2时，romberg返回了接近0的错误值，偏移下限后结果恢复正常。

问题分析 #

被积函数f(x)=x²cos²x是偶函数（满足f(-x)=f(x)），理论上∫_{-a}^{a}f(x)dx = 2∫_{0}^{a}f(x)dx，但Scipy 1.10.0的romberg算法在区间端点x=-pi/2处（此时cos(x)=0，函数值为0）出现了数值计算异常，导致迭代过程中出现错误的抵消或收敛判断，最终得到错误结果。

解决提示 #

• 利用偶函数对称性计算：直接通过两倍的半区间积分得到正确结果，这是最可靠的方式：

  res = 2 * romberg(f, 0, pi/2)
  print(res)  # 输出0.5065300316211875，与预期一致
  

• 轻微偏移区间端点：如果必须使用完整对称区间，可给端点添加极小的偏移量（如1e-8），避免触发算法异常：

  eps = 1e-8
  res = romberg(f, -pi/2 + eps, pi/2 - eps)
  

• 升级Scipy版本：该问题大概率是Scipy 1.10.0的版本bug，升级到1.11及以上版本后，romberg算法的相关问题可能已被修复。
• 替换为更稳定的积分方法：改用scipy.integrate.quad（自适应高斯-勒让德积分），该方法对这类积分场景的稳定性更好：

  from scipy.integrate import quad
  res, _ = quad(f, -pi/2, pi/2)
  print(res)  # 输出0.5065300316211875
  

内容的提问来源于stack exchange，提问作者Klaus Rohe