quadpy IntegrationError: Не удалось достичь допусков (abs: 1.49e-08, rel: 1.49e-08) с данными max_num_subintervals (= 50)

Question

Я использую 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)

Функция periodi c, это ее график:

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, Сципи все еще работает.

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.

Answer 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.

Answer 2

После еще нескольких экспериментов я обнаружил, что метод quadpy quad принимает тот же аргумент, что и метод scipys quad, который можно найти здесь:

• https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html#scipy .integrate.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

Однако в целом этот ответ отвечает на мой первоначальный вопрос.