Scipy dblquad计算结果异常，与Wolfram Alpha结果不符排查

问题重现 #

用户使用Scipy的dblquad计算二重积分，Python代码如下：

import numpy as np
from scipy.integrate import dblquad

def f(t,p):
  return np.sin(t) * np.exp( -760.45332604428756 *( -  (1.0 * 0.001) *  np.sin(t)*  np.cos(p) +   (-10 * 0.001/2) *  np.sin(t)**2 + (10 * 0.001/2)* np.sin(t)**2 * np.sin(p)**2 )   )

area = dblquad(f, 0, np.pi, -np.pi/2, np.pi/2)
print("area = ",area)

执行结果为：area = (6.355113632868125e-16, 2.1352999812415674e-07)，但用Wolfram Alpha计算同一积分：

Integrate [   Sin[t] Exp[ -  760.45332604428756  ( -  (1.0 * 0.001)  Sin[t] Cos[p ] +   (-10 * 0.001/2) Sin[t]^2 + (10 * 0.001/2) Sin[t]^2 Sin[p]^2 )  ], {t,0, Pi}, {p, -Pi/2, Pi/2}]

得到正确结果83.9147。

问题根源 #

Scipy的dblquad函数参数顺序与预期的积分顺序不匹配。

dblquad的定义规则是：dblquad(func, a, b, gfun, hfun)，其中：

• a和b是外层积分变量的上下限，对应func的第一个参数
• gfun和hfun是内层积分变量的上下限，对应func的第二个参数

你的代码中，dblquad(f, 0, np.pi, -np.pi/2, np.pi/2)等价于先对t（func第一个参数）在[0, π]积分，再对p（func第二个参数）在[-π/2, π/2]积分，但Wolfram Alpha的积分顺序是先对p积分，再对t积分，两者顺序完全相反。

更关键的是，原积分顺序下，内层对t积分时，指数部分会产生极大的正数值：
展开指数内的项：

-760.453 * [ -0.001 sin(t)cos(p) -0.005 sin²(t) +0.005 sin²(t) sin²(p) ]
= 760.453 * [0.001 sin(t)cos(p) +0.005 sin²(t)cos²(p)]

这个值是正数且量级很大，导致exp(-大正数)几乎为0，最终积分结果趋近于0。

修正方案 #

有两种方式可以修正：

方式1：调整dblquad的积分限顺序，匹配Wolfram的积分逻辑 #

把t的积分限设为内层，p的设为外层（内层积分限可以是外层变量的函数，这里是常数，直接用lambda表达式传入）：

area = dblquad(f, -np.pi/2, np.pi/2, lambda p: 0, lambda p: np.pi)

方式2：调换函数f的参数顺序，让内层积分变量对应func的第二个参数 #

修改函数参数为(p,t)，保持原积分限顺序不变：

def f(p,t):
  return np.sin(t) * np.exp( -760.45332604428756 *( -  (1.0 * 0.001) *  np.sin(t)*  np.cos(p) +   (-10 * 0.001/2) *  np.sin(t)**2 + (10 * 0.001/2)* np.sin(t)**2 * np.sin(p)**2 )   )

area = dblquad(f, 0, np.pi, -np.pi/2, np.pi/2)

两种方式执行后，都会得到接近83.9147的结果，与Wolfram Alpha的计算结果一致。

内容的提问来源于stack exchange，提问作者Ed Gan