0

我正在使用 quadpy 在 python 中集成一个函数。

功能

import numpy as np
T = 2*np.pi 

def ex1(t):
    return np.where(np.logical_and((t%T>=0), (t%T<np.pi)), t%T, np.pi)

该函数是周期性的，这是它的情节：

x = np.linspace(0, 6*T, 1000)
plt.plot(x, ex1(x))
plt.grid(True)
plt.show()

在此处输入图像描述

问题

我正在尝试集成此功能：

from scipy.integrate import quad
import quadpy

print(quadpy.quad(ex1, 0, 3))
print(quad(ex1, 0, 3))

生产

(array(4.5), array(1.41692995e-19))
(4.5, 4.9960036108132044e-14)

在 0 到 3 的时间间隔内，一切正常。但是，如果我将间隔增加到例如 4，scipy 仍然有效。

print(quad(ex1, 0, 4))

生产

(7.631568411183528, 1.0717732083155035e-08)

但

print(quadpy.quad(ex1, 0, 4))

生产

IntegrationError: Tolerances (abs: 1.49e-08, rel: 1.49e-08) could not be reached with the given max_num_subintervals (= 50).

问题

  • 如何防止此错误？我尝试添加一个名为的参数，max_num_subintervals但这似乎不起作用。
  • 我是否正确地使用了 quadpy 来做我想做的事情？我已经开始使用它，因为我想采用 scipy 不支持的复数的导​​数，并且我想有一个万能的解决方案，因此对于这些更简单的示例使用 quadpy scipy 就足够了。
4

2 回答 2

1

毕竟这是一个四倍的错误。现在修复：

import numpy as np
from scipy.integrate import quad
import quadpy

T = 2 * np.pi


def ex1(t):
    return np.where(np.logical_and((t % T >= 0), (t % T < np.pi)), t % T, np.pi)


print(quadpy.__version__)
print(quadpy.quad(ex1, 0, 4))
print(quad(ex1, 0, 4))

0.14.7
(array(7.63156841), array(1.13827355e-16))
(7.631568411183528, 1.0717732083161812e-08)

您可能使用的是过时的 quadpy 版本。

于 2020-04-14T14:21:00.397 回答
0

经过更多的实验，我发现 quadpy 的 quad 方法与 scipys quad 方法采用相同的参数，可以在这里找到：

使用epsabs, epsrel, limit可选参数，我可以防止错误：

print(quadpy.quad(ex1, 0, 4, epsabs=1e-1, epsrel=1e-1, limit=100))

生产

(array(7.6323447), array(0.01666253))

然而 scipys quad 方法的误差容限达到1-08，我真的无法用 quadpy 重现，即使将限制设置为非常高的值，例如1000000

这是为什么？

更新：看来我遇到的行为是一个错误，现在报告在：https ://github.com/nschloe/quadpy/issues/255

然而，总的来说，这个答案回答了我最初的问题。

于 2020-04-14T08:15:28.643 回答