scipy.integrate.

tanhsinh#

scipy.integrate.tanhsinh(f, a, b, *, args=(), log=False, maxlevel=None, minlevel=2, atol=None, rtol=None, preserve_shape=False, callback=None)[источник]#

Вычислить сходящийся интеграл численно с использованием квадратуры tanh-sinh.

На практике квадратура tanh-sinh достигает квадратичной сходимости для многих подынтегральных выражений: количество точных digits масштабируется примерно линейно с количеством вычислений функции [1].

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

Параметры:

fcallable

Функция для интегрирования. Сигнатура должна быть:

f(xi: ndarray, *argsi) -> ndarray

где каждый элемент xi является конечным вещественным числом и argsi является кортежем, который может содержать произвольное количество массивов, совместимых с xi. f должна быть поэлементной функцией: см. документацию параметра preserve_shape для деталей. Он не должен изменять массив xi или массивы в argsi. Если f возвращает значение с комплексным типом данных при вычислении в любой конечной точке, последующие аргументы x будет иметь комплексный тип данных (но с нулевой мнимой частью).

a, bfloat array_like

Действительные нижний и верхний пределы интегрирования. Должны быть совместимы для вещания друг с другом и с массивами в args. Элементы могут быть бесконечными.

argsкортеж из array_like, необязательно

Дополнительные позиционные аргументы в виде массивов для передачи в f. Массивы должны быть совместимы для трансляции друг с другом и массивами a и b. Если вызываемый объект, для которого требуется найти корень, требует аргументов, которые не транслируются с x, оберните этот вызываемый объект с f такой, что f принимает только x и поддерживающие вещание *args.

logbool, по умолчанию: False

Установка в True указывает, что f возвращает логарифм подынтегрального выражения и что atol и rtol выражаются как логарифмы абсолютной и относительной ошибок. В этом случае объект результата будет содержать логарифм интеграла и ошибки. Это полезно для подынтегральных функций, для которых численное исчезновение порядка или переполнение привело бы к неточностям. Когда log=True, подынтегральное выражение (экспонента от f) должен быть вещественным, но может быть отрицательным, в этом случае логарифм подынтегрального выражения является комплексным числом с мнимой частью, кратной нечетному числу π.

maxlevelint, по умолчанию: 10

Максимальный уровень уточнения алгоритма.

На нулевом уровне, f вызывается один раз, выполняя 16 вычислений функции. На каждом последующем уровне, f вызывается ещё раз, примерно удваивая количество выполненных вычислений функции. Соответственно, для многих подынтегральных выражений каждый последующий уровень будет удваивать количество точных цифр в результате (вплоть до ограничений точности с плавающей запятой).

Алгоритм завершится после завершения уровня maxlevel или после удовлетворения другого условия завершения, в зависимости от того, что наступит раньше.

minlevelint, по умолчанию: 2

Уровень, с которого начинается итерация (по умолчанию: 2). Это не изменяет общее количество вычислений функции или абсциссы, в которых функция вычисляется; это изменяет только количество раз f вызывается. Если minlevel=k, тогда подынтегральное выражение вычисляется на всех абсциссах от уровней 0 через k за один вызов. Обратите внимание, что если minlevel превышает maxlevel, предоставленный minlevel игнорируется, и minlevel устанавливается равным maxlevel.

atol, rtolfloat, опционально

Абсолютный допуск завершения (по умолчанию: 0) и относительный допуск завершения (по умолчанию: eps**0.75, где eps является точностью типа данных результата), соответственно. Должны быть неотрицательными и конечными, если log равно False и должно быть выражено как логарифм неотрицательного и конечного числа, если log равно True. Итерация остановится, когда res.error < atol или res.error < res.integral * rtol.

preserve_shapebool, по умолчанию: False

В следующем, "аргументы f” относится к массиву xi и любые массивы внутри argsi. Пусть shape будет трансляционной формой a, b, и все элементы args (что концептуально отличается от xi` and ``argsi передано в f).

Когда preserve_shape=False (по умолчанию), f должен принимать аргументы любой совместимые формы для вещания.
Когда preserve_shape=True, f должен принимать аргументы формы shape или shape + (n,), где (n,) это количество абсцисс, в которых вычисляется функция.

В любом случае, для каждого скалярного элемента xi[j] внутри xi, массив, возвращаемый f должен включать скаляр f(xi[j]) по тому же индексу. Следовательно, форма вывода всегда совпадает с формой ввода xi.

См. примеры.

callbackвызываемый объект, необязательный

Необязательная пользовательская функция, вызываемая перед первой итерацией и после каждой итерации. Вызывается как callback(res), где res является _RichResult аналогично тому, что возвращается tanhsinh (но содержащий текущие значения всех переменных на данной итерации). Если callback вызывает StopIteration, алгоритм немедленно завершится и tanhsinh вернет объект результата. callback не должен изменять res или его атрибуты.

Возвращает:

res_RichResult

Объект, похожий на экземпляр scipy.optimize.OptimizeResult со следующими атрибутами. (Описания написаны так, как будто значения будут скалярами; однако, если f возвращает массив, выходные данные будут массивами той же формы.)

successлогический массив

True когда алгоритм успешно завершился (статус 0). False в противном случае.

statusцелочисленный массив

Целое число, представляющее статус завершения алгоритма.

0 : Алгоритм сошелся к заданным допускам. -1 : (не используется) -2 : Достигнуто максимальное количество итераций. -3 : Встречено неконечное значение. -4 : Итерация была завершена callback. 1 : Алгоритм работает нормально (в callback только).

интегралмассив float

Оценка интеграла.

ошибкамассив float

Оценка ошибки. Доступно только если уровень два или выше был завершён; иначе NaN.

maxlevelцелочисленный массив

Максимальный уровень уточнения, используемый.

nfevцелочисленный массив

Количество точек, в которых f было вычислено.

Смотрите также

quad

Примечания

Реализует алгоритм, как описано в [1] с небольшими адаптациями для арифметики конечной точности, включая некоторые, описанные [2] и [3]. Схема tanh-sinh была первоначально представлена в [4].

Две схемы оценки ошибки описаны в [1] Раздел 5: одна попытка обнаружить и использовать квадратичную сходимость; другая просто сравнивает оценки интеграла на последовательных уровнях. Хотя ни одна из них не является теоретически строгой или консервативной, обе хорошо работают на практике. Наша оценка ошибки использует минимум этих двух схем с нижней границей eps * res.integral.

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

Ссылки

[1] (1,2,3)

Бейли, Дэвид Х., Картик Джеябалан и Сяое С. Ли. «Сравнение трёх схем квадратуры высокой точности». Experimental Mathematics 14.3 (2005): 317-329.

[2]

Vanherck, Joren, Bart Sorée, и Wim Magnus. «Tanh-sinh quadrature for single and multiple integration using floating-point arithmetic.» arXiv preprint arXiv:2007.15057 (2020).

[3]

van Engelen, Robert A. "Improving the Double Exponential Quadrature Tanh-Sinh, Sinh-Sinh and Exp-Sinh Formulas." https://www.genivia.com/files/qthsh.pdf

[4]

Takahasi, Hidetosi, and Masatake Mori. “Double exponential formulas for numerical integration.” Publications of the Research Institute for Mathematical Sciences 9.3 (1974): 721-741.

Примеры

Вычислить интеграл Гаусса:

>>> import numpy as np
>>> from scipy.integrate import tanhsinh
>>> def f(x):
...     return np.exp(-x**2)
>>> res = tanhsinh(f, -np.inf, np.inf)
>>> res.integral  # true value is np.sqrt(np.pi), 1.7724538509055159
1.7724538509055159
>>> res.error  # actual error is 0
4.0007963937534104e-16

Значение функции Гаусса (колоколообразной кривой) близко к нулю для аргументов, достаточно удалённых от нуля, поэтому значение интеграла по конечному интервалу почти одинаково.

>>> tanhsinh(f, -20, 20).integral
1.772453850905518

Однако при неблагоприятных пределах интегрирования схема интегрирования может не найти важную область.

>>> tanhsinh(f, -np.inf, 1000).integral
4.500490856616431

В таких случаях или при наличии особенностей внутри интервала разбейте интеграл на части с конечными точками в важных точках.

>>> tanhsinh(f, -np.inf, 0).integral + tanhsinh(f, 0, 1000).integral
1.772453850905404

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

>>> res = tanhsinh(f, 20, 30, rtol=1e-10)
>>> res.integral, res.error
(4.7819613911309014e-176, 4.670364401645202e-187)
>>> def log_f(x):
...     return -x**2
>>> res = tanhsinh(log_f, 20, 30, log=True, rtol=np.log(1e-10))
>>> np.exp(res.integral), np.exp(res.error)
(4.7819613911306924e-176, 4.670364401645093e-187)

Пределы интегрирования и элементы args могут быть транслируемыми массивами, и интегрирование выполняется поэлементно.

>>> from scipy import stats
>>> dist = stats.gausshyper(13.8, 3.12, 2.51, 5.18)
>>> a, b = dist.support()
>>> x = np.linspace(a, b, 100)
>>> res = tanhsinh(dist.pdf, a, x)
>>> ref = dist.cdf(x)
>>> np.allclose(res.integral, ref)
True

По умолчанию, preserve_shape равно False, и поэтому вызываемый объект f может вызываться с массивами любых совместимых форм. Например:

>>> shapes = []
>>> def f(x, c):
...    shape = np.broadcast_shapes(x.shape, c.shape)
...    shapes.append(shape)
...    return np.sin(c*x)
>>>
>>> c = [1, 10, 30, 100]
>>> res = tanhsinh(f, 0, 1, args=(c,), minlevel=1)
>>> shapes
[(4,), (4, 34), (4, 32), (3, 64), (2, 128), (1, 256)]

Чтобы понять, откуда берутся эти формы - и лучше понять, как tanhsinh вычисляет точные результаты - обратите внимание, что более высокие значения c соответствуют синусоидам более высокой частоты. Синусоиды более высокой частоты усложняют подынтегральное выражение, поэтому для достижения целевой точности требуется больше вычислений функции:

>>> res.nfev
array([ 67, 131, 259, 515], dtype=int32)

Начальный shape, (4,), соответствует вычислению интеграла в одной абсциссе и всех четырех частотах; это используется для проверки входных данных и определения размера и типа данных массивов, которые хранят результаты. Следующая форма соответствует вычислению интеграла в начальной сетке абсцисс и всех четырех частотах. Последующие вызовы функции удваивают общее количество абсцисс, в которых функция была вычислена. Однако в последующих вычислениях функции интеграл вычисляется на меньшем количестве частот, потому что соответствующий интеграл уже сошелся к требуемой точности. Это экономит вычисления функции для повышения производительности, но требует, чтобы функция принимала аргументы любой формы.

«Векторно-значные» подынтегральные выражения, такие как те, что используются с scipy.integrate.quad_vec, вряд ли удовлетворяют этому требованию. Например, рассмотрим

>>> def f(x):
...    return [x, np.sin(10*x), np.cos(30*x), x*np.sin(100*x)**2]

Этот подынтегральный выражение несовместим с tanhsinh как написано; например, форма вывода не будет такой же, как форма x. Такая функция мог может быть преобразовано в совместимую форму с введением дополнительных параметров, но это было бы неудобно. В таких случаях более простое решение — использовать preserve_shape.

>>> shapes = []
>>> def f(x):
...     shapes.append(x.shape)
...     x0, x1, x2, x3 = x
...     return [x0, np.sin(10*x1), np.cos(30*x2), x3*np.sin(100*x3)]
>>>
>>> a = np.zeros(4)
>>> res = tanhsinh(f, a, 1, preserve_shape=True)
>>> shapes
[(4,), (4, 66), (4, 64), (4, 128), (4, 256)]

Здесь трансляционная форма a и b является (4,). С preserve_shape=True, функция может быть вызвана с аргументом x формы (4,) или (4, n), и это то, что мы наблюдаем.