scipy.integrate.

кубатура#

scipy.integrate.кубатура(f, a, b, *, правило='gk21', rtol=1e-08, atol=0, максимальное_количество_подразделений=10000, args=(), workers=1, точки=None)[источник]#

Адаптивное кубатурное интегрирование многомерной функции со значениями в виде массивов.

Для произвольного правила интегрирования эта функция возвращает оценку интеграла с требуемой точностью по области, заданной массивами a и b указывая углы гиперкуба.

Сходимость не гарантируется для всех интегралов.

Параметры:

fcallable

Функция для интегрирования. f должен иметь сигнатуру:

f(x : ndarray, *args) -> ndarray

f должен принимать массивы x формы:

(npoints, ndim)

и выходные массивы формы:

(npoints, output_dim_1, ..., output_dim_n)

В этом случае, cubature вернет массивы формы:

(output_dim_1, ..., output_dim_n)

a, barray_like

Нижний и верхний пределы интегрирования как одномерные массивы, задающие левые и правые концы интервалов, по которым происходит интегрирование. Пределы могут быть бесконечными.

правилоstr, optional

Правило, используемое для оценки интеграла. Если передается строка, варианты: "gauss-kronrod" (21 узел) или "genz-malik" (степень 7). Если правило типа "gauss-kronrod" указано для n-мерный интегранд, используется соответствующее правило декартова произведения. “gk21”, “gk15” также поддерживаются для совместимости с quad_vec. См. Примечания.

rtol, atolfloat, опционально

Относительные и абсолютные допуски. Итерации выполняются до тех пор, пока ошибка не будет оценена как меньшая, чем atol + rtol * abs(est). Здесь rtol управляет относительной точностью (количество верных цифр), в то время как atol управляет абсолютной точностью (количество правильных десятичных знаков). Для достижения желаемой rtol, установите atol быть меньше наименьшего значения, которое можно ожидать от rtol * abs(y) так что rtol преобладает над допустимой погрешностью. Если atol больше, чем rtol * abs(y) количество верных цифр не гарантируется. И наоборот, для достижения желаемого atol, установите rtol такой, что rtol * abs(y) всегда меньше, чем atol. Значения по умолчанию: 1e-8 для rtol и 0 для atol.

максимальное_количество_подразделенийint, необязательный

Верхняя граница количества подразделений для выполнения. По умолчанию 10,000.

argsкортеж, необязательный

Дополнительные позиционные аргументы, передаваемые в f, если есть.

workersint или вызываемый объект, подобный отображению, необязательный

Если workers является целым числом, часть вычислений выполняется параллельно, разделенная на это количество задач (используя multiprocessing.pool.Pool). Укажите -1 использовать все доступные ядра процессора. В качестве альтернативы можно передать вызываемый объект, подобный map, например multiprocessing.pool.Pool.map для параллельной оценки популяции. Эта оценка выполняется как workers(func, iterable).

точкисписок array_like, опционально

Список точек, которые следует избегать при вычислении f при условии, что используемое правило не вычисляет f на границе области (что имеет место для всех правил Генца-Малика и Гаусса-Кронрода). Это может быть полезно, если f имеет особенность в указанной точке. Это должен быть список массивоподобных объектов, где каждый элемент имеет длину ndim. По умолчанию пусто. См. Примеры.

Возвращает:

resobject

Объект, содержащий результаты оценки. Имеет следующие атрибуты:

оценкаndarray: Оценка значения интеграла по всей указанной области.
ошибкаndarray: Оценка ошибки аппроксимации по всей указанной области.
statusstr: Успешна ли оценка. Может быть: "converged", "not_converged".
подразделенияint: Количество выполненных подразделений.
atol, rtolfloat: Запрошенные допуски для аппроксимации.
регионы: список объектов: Список объектов, содержащих оценки интеграла по меньшим областям определения.

Каждый объект в regions имеет следующие атрибуты:

a, bndarray: Точки, описывающие углы области. Если исходный интеграл содержал бесконечные пределы или был по области, описанной область, затем a и b находятся в преобразованных координатах.
оценкаndarray: Оценка значения интеграла по этой области.
ошибкаndarray: Оценка ошибки аппроксимации в этой области.

Примечания

Алгоритм использует аналогичный алгоритм quad_vec, который сам основан на реализации алгоритмов DQAG* из QUADPACK, реализующей глобальный контроль ошибок и адаптивное разбиение.

Источник узлов и весов, используемых для квадратуры Гаусса-Кронрода, можно найти в [1], и алгоритм вычисления узлов и весов в кубатуре Генца-Малика можно найти в [2].

Правила, в настоящее время поддерживаемые через правило аргумента:

"gauss-kronrod", 21-узловая Гаусса-Кронрода
"genz-malik", n-узловая Genz-Malik

Если использовать Гаусса-Кронрода для параметра, результат будет корректным только если он соответствует n-мерный интегранд, где n > 2, тогда соответствующее правило декартова произведения будет найдено путем взятия декартова произведения узлов в одномерном случае. Это означает, что количество узлов масштабируется экспоненциально как 21^n в случае Гаусса-Кронрода, что может быть проблематичным в умеренном числе измерений.

Genz-Malik обычно менее точен, чем Gauss-Kronrod, но имеет гораздо меньше узлов, поэтому в этой ситуации использование "genz-malik" может быть предпочтительнее.

Бесконечные пределы обрабатываются с помощью соответствующего преобразования переменных. Предполагая a = [a_1, ..., a_n] и b = [b_1, ..., b_n]:

Если \(a_i = -\infty\) и \(b_i = \infty\), i-я переменная интегрирования будет использовать преобразование \(x = \frac{1-|t|}{t}\) и \(t \in (-1, 1)\).

Если \(a_i \ne \pm\infty\) и \(b_i = \infty\), i-я переменная интегрирования будет использовать преобразование \(x = a_i + \frac{1-t}{t}\) и \(t \in (0, 1)\).

Если \(a_i = -\infty\) и \(b_i \ne \pm\infty\), i-я переменная интегрирования будет использовать преобразование \(x = b_i - \frac{1-t}{t}\) и \(t \in (0, 1)\).

Ссылки

[1]

R. Piessens, E. de Doncker, Quadpack: A Subroutine Package for Automatic Integration, files: dqk21.f, dqk15.f (1983).

[2]

A.C. Genz, A.A. Malik, Remarks on algorithm 006: An adaptive algorithm for numerical integration over an N-dimensional rectangular region, Journal of Computational and Applied Mathematics, Volume 6, Issue 4, 1980, Pages 295-302, ISSN 0377-0427 DOI:10.1016/0771-050X(80)90039-X

Примеры

1D интеграл с векторным выводом:

\[\int^1_0 \mathbf f(x) \text dx\]

Где f(x) = x^n и n = np.arange(10) является вектором. Поскольку правило не указано, используется значение по умолчанию «gk21», что соответствует интегрированию Гаусса-Кронрода с 21 узлом.

>>> import numpy as np
>>> from scipy.integrate import cubature
>>> def f(x, n):
...    # Make sure x and n are broadcastable
...    return x[:, np.newaxis]**n[np.newaxis, :]
>>> res = cubature(
...     f,
...     a=[0],
...     b=[1],
...     args=(np.arange(10),),
... )
>>> res.estimate
 array([1.        , 0.5       , 0.33333333, 0.25      , 0.2       ,
        0.16666667, 0.14285714, 0.125     , 0.11111111, 0.1       ])

7D интеграл с выходным массивом произвольной формы:

f(x) = cos(2*pi*r + alphas @ x)

для некоторого r и alphas, и интеграл вычисляется по единичному гиперкубу, \([0, 1]^7\). Поскольку интеграл находится в умеренном количестве измерений, используется 'genz-malik', а не 'gauss-kronrod' по умолчанию, чтобы избежать построения правила произведения с \(21^7 \approx 2 \times 10^9\) узлы.

>>> import numpy as np
>>> from scipy.integrate import cubature
>>> def f(x, r, alphas):
...     # f(x) = cos(2*pi*r + alphas @ x)
...     # Need to allow r and alphas to be arbitrary shape
...     npoints, ndim = x.shape[0], x.shape[-1]
...     alphas = alphas[np.newaxis, ...]
...     x = x.reshape(npoints, *([1]*(len(alphas.shape) - 1)), ndim)
...     return np.cos(2*np.pi*r + np.sum(alphas * x, axis=-1))
>>> rng = np.random.default_rng()
>>> r, alphas = rng.random((2, 3)), rng.random((2, 3, 7))
>>> res = cubature(
...     f=f,
...     a=np.array([0, 0, 0, 0, 0, 0, 0]),
...     b=np.array([1, 1, 1, 1, 1, 1, 1]),
...     rtol=1e-5,
...     rule="genz-malik",
...     args=(r, alphas),
... )
>>> res.estimate
 array([[-0.79812452,  0.35246913, -0.52273628],
        [ 0.88392779,  0.59139899,  0.41895111]])

Параллельные вычисления с workers:

>>> from concurrent.futures import ThreadPoolExecutor
>>> with ThreadPoolExecutor() as executor:
...     res = cubature(
...         f=f,
...         a=np.array([0, 0, 0, 0, 0, 0, 0]),
...         b=np.array([1, 1, 1, 1, 1, 1, 1]),
...         rtol=1e-5,
...         rule="genz-malik",
...         args=(r, alphas),
...         workers=executor.map,
...      )
>>> res.estimate
 array([[-0.79812452,  0.35246913, -0.52273628],
        [ 0.88392779,  0.59139899,  0.41895111]])

2D интеграл с бесконечными пределами:

\[\int^{ \infty }_{ -\infty } \int^{ \infty }_{ -\infty } e^{-x^2-y^2} \text dy \text dx\]

>>> def gaussian(x):
...     return np.exp(-np.sum(x**2, axis=-1))
>>> res = cubature(gaussian, [-np.inf, -np.inf], [np.inf, np.inf])
>>> res.estimate
 3.1415926

1D интеграл с избеганием сингулярностей с использованием точки:

\[\int^{ 1 }_{ -1 } \frac{\sin(x)}{x} \text dx\]

Необходимо использовать точки параметр для избежания вычисления f в начале координат.

>>> def sinc(x):
...     return np.sin(x)/x
>>> res = cubature(sinc, [-1], [1], points=[[0]])
>>> res.estimate
 1.8921661