使用Sympy和Scipy的dblquad导致数值积分无法收敛

Question

简介

Python Scipy 库提供了

dblquad

函数，可用于数值计算二维积分。 Python Sympy 库提供了在 Python 中定义和计算符号表达式的工具。

将

dblquad

函数与 Sympy 表达式组合时，会发生如下所述的不良行为。

我的库版本和设置

• Python 3.9.13，
• Scipy 1.13.1
• 交感1.12.1.

下面的所有代码都假定以导入为前缀

import numpy as np

from scipy.integrate import dblquad
from sympy.functions.special.beta_functions import beta

问题陈述

运行以下代码时：

function = lambda x, y: beta(2, 2) if 0 < x < 1 and 0 < y < 1 else 0

result = dblquad(lambda x, y: function(y, x),
    0.99 / 2, 1,  # boundaries for x
    lambda x: 0.99 / x, lambda x: np.inf  # boundaries for y
)

程序不会终止（至少在我运行代码的合理时间内不会终止）。这让我感到惊讶，原因如下。有人解释一下吗

将 Sympy 的 beta(2, 2) 替换为数值 1/6 时，程序几乎立即运行

如果我们运行的不是上面的代码

function = lambda x, y: 1 / 6 if 0 < x < 1 and 0 < y < 1 else 0

result = dblquad(lambda x, y: function(y, x),
    0.99 / 2, 1,  # boundaries for x
    lambda x: 0.99 / x, lambda x: np.inf  # boundaries for y
)

（即我们只需将 beta(2, 2) 替换为数值 1/6），程序很快就会终止并得到结果。

在整个第一象限上积分时，程序有效

如果我们运行的不是初始代码

function = lambda x, y: beta(2, 2) if 0 < x < 1 and 0 < y < 1 else 0

result = dblquad(lambda x, y: function(y, x),
    0, np.inf,  # boundaries for x
    0, np.inf  # boundaries for y
)

代码运行速度非常快。

最后一个问题

为什么会出现上述行为以及如何让初始程序运行。 （在我想到的应用程序中，我需要 Sympy 的更复杂的行为，我无法像用 1/6 替换 beta(2, 2) 那样轻松地用其他东西替换它）。

Answer 1

虽然

dblquad

的文档没有明确说明被积数的返回类型，但假设它必须是数字并不容易。

type(beta(2, 2))  # beta

看起来不是数字。

import numpy as np
np.issubdtype(float, np.number)  # True
np.issubdtype(np.float64, np.number)  # True
np.issubdtype(type(beta(2, 2)), np.number)  # False

不是；无论如何，不是根据 NumPy 对数字的定义。而且 SciPy 与 NumPy 有着密切的联系，因此 NumPy 对数字的定义是一个合理的使用方式。

所以我建议，代码曾经有效比它并不总是有效更令人惊讶。

现在，如果它曾经确实起作用，那可能是因为在某个时候这个

beta

对象被强制转换为浮点数。例如，类似：

float(beta(2, 2))  # 0.16666666666666666

在某个时刻发生。

将 Sympy 的 beta(2, 2) 替换为其数值 1/6 时，程序几乎立即运行

使用 SymPy 评估

beta(2, 2)

然后转换为

float

比评估

1/6

慢得多。

%timeit 1/6
# 15.3 ns ± 3.47 ns per loop (mean ± std. dev. of 7 runs, 100000000 loops each)
%timeit float(beta(2, 2))
# 163 µs ± 31.1 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

因此，如果积分确实运行的话，可以合理地预期积分运行时间将至少延长 10000 倍。

在整个第一象限上积分时，该程序有效

在 Colab 上，两者都有效，但整个第一象限是一个简单得多的问题（具有不同的答案），因此执行速度要快得多。

%timeit dblquad(lambda x, y: function(y, x), 0, np.inf, 0, np.inf)
# 78.7 ms ± 15.1 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit dblquad(lambda x, y: function(y, x), 0.99 / 2, 1, lambda x: 0.99 / x, lambda x: np.inf)
# 9.92 s ± 438 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

即使您使用

1/6

而不是

beta(2, 2)

，在整个第一象限上积分也会快得多。

如何让初始程序运行

有一个关于

lambdify

的很好的评论：

import numpy as np

from scipy.integrate import dblquad
from sympy.functions.special.beta_functions import beta

from sympy import lambdify
from sympy.abc import x, y, z

# first argument: a sequence of symbolic arguments
# second argument: symbolic expression
# output: function that accepts numbers and returns numbers
my_beta = lambdify((x, y), beta(2, 2))  # the expression `beta(2, 2)` doesn't use `x` or `y`, but it could be replaced with an expression that does

现在我们可以将数字参数传递给

my_beta

，它将返回

beta(2, 2)

作为数字：

my_beta(1, 10)  # 0.16666666666666666
%timeit my_beta(1, 10)
# 1.98 µs ± 64.6 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

因此您可以在代码中使用它，例如：

function = lambda x, y: my_beta(x, y) if 0 < x < 1 and 0 < y < 1 else 0

result = dblquad(lambda x, y: function(y, x),
    0.99 / 2, 1,  # boundaries for x
    lambda x: 0.99 / x, lambda x: np.inf  # boundaries for y
)

另一种选择是使用

scipy

特殊函数直接编写表达式。例如，使用

scipy.special.beta

来评估 beta 函数。

from scipy.special import beta
function = lambda x, y: beta(2, 2) if 0 < x < 1 and 0 < y < 1 else 0
result = dblquad(lambda x, y: function(y, x),
    0.99 / 2, 1,  # boundaries for x
    lambda x: 0.99 / x, lambda x: np.inf  # boundaries for y
)