Как quad в SciPy интегрирует бесконечные интервалы: qagie, отображение на [0,1] и практические приёмы

Объясняем, как SciPy quad интегрирует бесконечные пределы через qagie из QUADPACK: отображение на [0,1], вес T^-2, симметрия для (-∞,∞), работа с points.

Интегрирование функций на бесконечных интервалах кажется обманчиво простым, когда используешь quad в SciPy, но под капотом работает точная механика. Понимание координатного отображения, которое применяет процедура qagie из QUADPACK, помогает рассуждать о численной устойчивости, сходимости и о том, как работать с такими возможностями, как points, на конечной области.

Постановка задачи

Ниже — минимальный пример, который интегрирует стандартную нормальную плотность (PDF) от 0 до бесконечности с помощью SciPy. Цель — понять, как quad обрабатывает бесконечные пределы.

import scipy
import numpy as np
def pdf_normal(z):
    return np.exp(-(z**2) / 2) / np.sqrt(2 * np.pi)
val, err, dbg = scipy.integrate.quad(pdf_normal, 0.0, np.inf, full_output=True)
print(val)

В документации SciPy сказано, что бесконечные пределы обрабатываются с помощью преобразования координат через qagie из QUADPACK. Их формулировка:

qagie

выполняет интегрирование на бесконечных интервалах. Бесконечный диапазон отображается на конечный интервал, после чего применяется та же стратегия, что и в QAGS.

Интерес представляет точное отображение на конечный интервал [0, 1] и то, что оно означает для формы подынтегральной функции и численного поведения.

Что на самом деле происходит под капотом

Отображение, используемое в qagie, задано явно. Для полу-бесконечного интеграла с нижним пределом A выполняется подстановка T из (0, 1], отображаемая в X из [A, ∞) по формуле X = A + (1 − T)/T. Это преобразование сводит исходный интеграл к интегралу по единичному отрезку с модифицированной подынтегральной функцией. Для двух сторон бесконечности та же идея применяется симметрично: для каждого T вычисляют функцию как в точке X, так и в точке −X.

Подробности описаны в: Piessens, Robert; de Doncker-Kapenga, Elise; Überhuber, Christoph W.; Kahaner, David (1983). QUADPACK: A subroutine package for automatic integration. Springer-Verlag. Стр. 64–65.

Один бесконечный предел: [A, ∞)

Используется подстановка X = A + (1 − T)/T при T в (0, 1]. Когда T → 1, X → A. Когда T → 0+, X → ∞. Дифференцирование даёт dX/dT = −T^(−2). Переворот пределов убирает минус, в результате получаем интеграл по T от 0 до 1 с подынтегральной функцией F(A + (1 − T)/T) · T^(−2). Это суть того, как qagie переносит полу-бесконечный интеграл на [0, 1].

Два бесконечных предела: (−∞, ∞)

Для интегралов на всей прямой применяется та же схема, но обе «хвостовые» части обрабатываются симметрично. Каждому T из (0, 1] соответствует X = (1 − T)/T, и подынтегральная функция берётся как F(X) и F(−X), складываясь в F(X) + F(−X), снова с весом T^(−2).

Практическое решение: ручная репараметризация

Вы можете выполнить ту же подстановку сами, интегрировать по [0, 1] и посмотреть, как ведёт себя преобразованная подынтегральная функция. Это полезно для диагностики и для использования опции points, которая иначе несовместима с бесконечными пределами.

import scipy
import numpy as np
import matplotlib.pyplot as plt
def pdf_normal(z):
    return np.exp(-(z**2) / 2) / np.sqrt(2 * np.pi)
def reparam_semi_infinite(fn, a):
    """Map T in (0, 1] to X in [a, +inf) via X = a + (1 - T)/T."""
    def wrapped(t):
        return fn(a + (1 - t) / t) * (t ** -2)
    return wrapped
def reparam_full_line(fn):
    """Map T in (0, 1] to X in (-inf, +inf) using symmetry."""
    def wrapped(t):
        x = (1 - t) / t
        return (fn(x) + fn(-x)) * (t ** -2)
    return wrapped
# Полу-бесконечный случай: интегрируем нормальную PDF от 0 до +inf через T в [eps, 1]
pdf_mapped = reparam_semi_infinite(pdf_normal, a=0.0)
# Преобразованная подынтегральная функция имеет особенность при T=0 (соответствует X = +inf),
# поэтому стартуем чуть правее нуля.
t0 = 1e-8
res, err, meta = scipy.integrate.quad(pdf_mapped, t0, 1.0, full_output=True)
print(res)
# Визуализируем исходную и преобразованную подынтегральные функции
xs1 = np.linspace(0.0, 5.0, 100)
plt.title("Normal PDF in X domain")
plt.plot(xs1, pdf_normal(xs1))
plt.show()
xs2 = np.linspace(t0, 1.0, 100)
plt.title("Transformed integrand in T domain")
plt.plot(xs2, pdf_mapped(xs2))
plt.show()

Интеграл по T не определён при T = 0, что под подстановкой соответствует X = ∞; поэтому нижнюю границу для T нужно брать как малое положительное число. В этой схеме вы интегрируете по [t0, 1] и получаете точно то же значение, что и исходный полу-бесконечный интеграл.

Если нужно задать разрывы/точки излома подынтегральной функции, теперь можно использовать аргумент points, потому что область — [0, 1]. Точки, заданные в исходной X-области, переносятся в T по формуле t = 1 / (1 − a + x), где a — нижний предел исходного интеграла.

Зачем это нужно

Знание отображения проясняет три практических момента. Во-первых, сингулярность при T = 0 — артефакт преобразования и соответствует поведению на бесконечности в исходной переменной; это помогает понять, почему интегратор испытывает трудности или почему нельзя брать T = 0 точно. Во-вторых, осмотр преобразованной подынтегральной функции показывает, не приходится ли интегратору бороться с резкими изменениями возле концов отрезка. В-третьих, как только интеграл живёт на [0, 1], становятся доступны такие возможности, как points, что может повысить надёжность возле резких изменений после переноса точек излома по данной формуле.

Вывод

quad обрабатывает бесконечные диапазоны, перепараметризуя их на [0, 1]. Для [A, ∞) используется X = A + (1 − T)/T с весом T^(−2); для (−∞, ∞) обе «стороны» берутся симметрично как F(X) + F(−X) с тем же весом. Вы можете воспроизвести преобразование вручную, чтобы визуализировать подынтегральную функцию, расставить точки излома по формуле t = 1/(1 − a + x) и лучше понимать, где возникают численные сложности. Если проблемы наблюдаются возле концов по T, помните, что это отражает поведение на бесконечности в исходной переменной, и соответственно подберите пределы интегрирования или диагностические приёмы.