Стабилизация интегралов в модели Томаса—Хопфилда: points в quad, точность и хвост к бесконечности

Разбираем, почему scipy.integrate.quad чувствителен к r_max в модели Томаса—Хопфилда, и как стабилизировать интеграл: points, epsabs=0, epsrel и разбиение хвоста к бесконечности.

Когда физическую модель с двумя несобственными интегралами переносишь в код, мелочи численного характера способны как улучшить, так и испортить график. Именно это и происходит с моделью Томаса—Хопфилда при вычислении через scipy.integrate.quad: результат чересчур зависит от верхней границы интегрирования, а векторизация для построения графика лишь усиливает непонимание. Хорошая новость в том, что такое поведение напрямую следует из устройства адаптивной квадратуры и его можно корректно исправить.

Как воспроизвести проблему

Модель интегрирует два ядра по r от 0 до бесконечности и собирает из них зависящий от времени сигнал. Прямое интегрирование до бесконечности неудобно, поэтому используют конечную отсечку r_max. Чувствительность к r_max — это симптом, на который и стоит смотреть.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

# Скорость рекомбинации (экспоненциально убывает с расстоянием)
def rate_fn(r, w0, r0):
    return w0 * np.exp((-2 * r) / r0)

# Первое ядро модели

def kern_a(r, w0, r0, t):
    return rate_fn(r, w0, r0) * np.exp(-rate_fn(r, w0, r0) * t) * r**2

# Второе ядро модели (используется в показателе экспоненты)

def kern_b(r, w0, r0, t):
    return (np.exp(-rate_fn(r, w0, r0) * t) - 1) * r**2

# Модель Томаса—Хопфилда, составленная из двух интегралов

def hopfield_fn(t, w0, r0, dens, amp, bias):
    r_cap = 250e-7  # прагматическая отсечка вместо бесконечности
    val1 = quad(kern_a, 0, r_cap, args=(w0, r0, t))[0]
    val2 = quad(kern_b, 0, r_cap, args=(w0, r0, t))[0]
    return amp * ((4 * np.pi * dens * val1) * np.exp(4 * np.pi * dens * val2)) + bias

# Векторизованная обёртка, принимающая массив t
hopfield_vec = np.vectorize(hopfield_fn)

# Пример набора параметров
w0 = 1e6
r0 = 1.2e-9
dens = 3e-19
amp = 1
bias = 0

# Теоретическое значение при t=0 по формуле из литературы
teor_t0 = dens * np.pi * w0 * r0**3
calc_t0 = hopfield_fn(0, w0, r0, dens, amp, bias)
print('Theoretical value for t=0:', teor_t0)
print('Calculated value for t=0:', calc_t0)

# График на логарифмических осях по x и y
T = np.linspace(0, 100000000, 100000)
Y = hopfield_vec(T, w0, r0, dens, amp, bias)
plt.plot(T, Y, color='r', linestyle='-', linewidth=0.7)
plt.yscale('log')
plt.xscale('log')
plt.show()

При изменении r_cap кривая растягивается, и совпадение при t=0 порой лучше для меньшего r_cap, чем для большего. Это не случайность: так и должен вести себя адаптивный интегратор, если ему не указать, где у под integrand есть сложный участок.

Как на самом деле работает quad и почему это важно

Пользуясь scipy.integrate.quad впервые, я задаюсь вопросом, как quad вычисляет результат (в каком порядке он проходит по массиву, считает интегралы и возвращает значения).

quad применяет адаптивную схему Гаусса—Кронрода. Он берет значения подынтегральной функции в наборе узлов, строит на текущем подынтервале две оценки — низкого и высокого порядка, — сравнивает их, чтобы оценить локальную ошибку, и уточняет сетку только там, где расхождение намекает на «интересное» поведение. Если ни один из узлов не попадает близко к пику или области резких изменений, алгоритм может ошибочно решить, что функция там почти плоская, и пойти дальше.

Минимальный пример делает это наглядным. Интегрирование стандартной нормальной плотности, центрированной в 0, по широкой симметричной области дает 1 — как и ожидается.

import scipy
import numpy as np

pdf0 = scipy.stats.norm(loc=0, scale=1)
val, err = scipy.integrate.quad(pdf0.pdf, -1000, 1000)
print('result', val, 'error', err)

Теперь сдвинем ту же распределение так, чтобы центр оказался в 100, и проинтегрируем на тех же пределах. Теоретически ничего не меняется. На практике, если адаптивный семплер ни разу не заглянет в окрестность x=100, он решит, что функция везде ничтожна, и вернет ноль.

pdf100 = scipy.stats.norm(loc=100, scale=1)
val, err = scipy.integrate.quad(pdf100.pdf, -1000, 1000)
print('result', val, 'error', err)

Чтобы это устранить, нужно подсказать quad, где у функции есть структура. Аргумент points делает ровно это — принудительно разбивает область вокруг указанных координат.

val, err = scipy.integrate.quad(pdf100.pdf, -1000, 1000, points=[100])
print('result', val, 'error', err)

Возвращаясь к модели Томаса—Хопфилда: у интегрируемых выражений есть чувствительные к масштабу особенности при малых r — примерно в районе от 1e-9 до 1e-7 в приведенном примере. Если не подтолкнуть интегратор смотреть именно туда, результаты сильно зависят от произвольно выбранных границ. Есть и второй момент. Стандартная абсолютная точность у quad рассчитана на величины порядка единицы. При крайне малых значениях типичный epsabs оказывается слишком велик. Установка epsabs в ноль и переход к относительной точности задают интегратору адекватные цели по погрешности.

Надёжное решение для конечной границы

В этом варианте мы явно указываем quad, где «происходит действие», и ужесточаем допуски. Единственное изменение — более точные и стабильные результаты.

def hopfield_fn(t, w0, r0, dens, amp, bias):
    r_cap = 1.0
    anchors = np.array([1e-9, 1e-8, 1e-7])
    val1 = quad(kern_a, 0, r_cap, args=(w0, r0, t), points=anchors, epsabs=0, epsrel=1.49e-8)[0]
    val2 = quad(kern_b, 0, r_cap, args=(w0, r0, t), points=anchors, epsabs=0, epsrel=1.49e-8)[0]
    return amp * ((4 * np.pi * dens * val1) * np.exp(4 * np.pi * dens * val2)) + bias

Это сразу убирает чувствительность к r_cap внутри конечного интервала и согласует значение при t=0 с теоретической формулой — при условии, что интеграл на интересующем промежутке посчитан достаточно точно.

Расширение на интегралы по [0, ∞)

Аргумент points несовместим с бесконечными пределами в используемой библиотеке. Обойти это просто: разделите интеграл на конечную часть, где используете points, и «хвост», уходящий к бесконечности, без них. На практике выбирайте отсечку так, чтобы хвост был пренебрежимо мал или хотя бы сглажен.

def hopfield_fn(t, w0, r0, dens, amp, bias):
    r_cap = 1.0
    anchors = np.array([1e-9, 1e-8, 1e-7])
    a_finite = quad(kern_a, 0, r_cap, args=(w0, r0, t), points=anchors, epsabs=0, epsrel=1.49e-8)[0]
    a_tail   = quad(kern_a, r_cap, np.inf, args=(w0, r0, t), epsabs=0, epsrel=1.49e-8)[0]
    b_finite = quad(kern_b, 0, r_cap, args=(w0, r0, t), points=anchors, epsabs=0, epsrel=1.49e-8)[0]
    b_tail   = quad(kern_b, r_cap, np.inf, args=(w0, r0, t), epsabs=0, epsrel=1.49e-8)[0]
    a_val = a_finite + a_tail
    b_val = b_finite + b_tail
    return amp * ((4 * np.pi * dens * a_val) * np.exp(4 * np.pi * dens * b_val)) + bias

С этим изменением теоретическая проверка при t=0 сходится на [0, ∞) следующим образом.

Theoretical value for t=0: 1.628601631620949e-39
Calculated value for t=0: 1.628601631620949e-39

О np.vectorize и альтернативах

np.vectorize не является источником интеграционной проблемы. Это удобная обёртка, которая по сути просто итерируется по t и собирает результаты. Прямой цикл на Python эквивалентен и порой проще для отладки.

grid_t = np.linspace(0, 1e8, 1000)
vals = np.array([hopfield_fn(tt, w0, r0, dens, amp, bias) for tt in grid_t])

Сам quad — скалярный интегратор: его подынтегральная функция вычисляется при скалярных r, так что передавать t по одному — правильный способ работы.

Опционально: как автоматически находить «интересные» r

Когда заранее сложно угадать хорошие опорные точки, один из программных подходов — искать нетривиальные стационарные точки подынтегральных функций и использовать их для разбиения области. Ниже программа символически дифференцирует каждое ядро, находит корень первой производной вдали от r=0 и по нему расставляет опоры. Затем эти опоры передаются в quad с ужатыми допусками. Это решение сложнее, его приводим как автоматизированную альтернативу.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
import sympy
import mpmath
import scipy
import math
from collections import defaultdict
mpmath.mp.dps = 50

def rate_fn(r, w0, r0):
    return w0 * np.exp((-2 * r) / r0)

def kern_a(r, w0, r0, t):
    return rate_fn(r, w0, r0) * np.exp(-rate_fn(r, w0, r0) * t) * r**2

def kern_b(r, w0, r0, t):
    return (np.exp(-rate_fn(r, w0, r0) * t) - 1) * r**2

w0_symb = sympy.Symbol('w0')
r0_symb = sympy.Symbol('r0')
t_symb = sympy.Symbol('t')
r_symb = sympy.Symbol('r')
w0c = 1e6
r0c = 1.2e-9

def rate_expr():
    return w0_symb * sympy.exp((-2 * r_symb) / r0_symb)

def kern_a_expr():
    return rate_expr() * sympy.exp(-rate_expr() * t_symb) * r_symb**2

def kern_b_expr():
    return (sympy.exp(-rate_expr() * t_symb) - 1) * r_symb**2

def expr_derivative(expr, w0, r0, n=1):
    expr = expr.subs(w0_symb, w0)
    expr = expr.subs(r0_symb, r0)
    deriv = sympy.diff(expr, r_symb, n)
    deriv = sympy.simplify(deriv)
    return sympy.lambdify((r_symb, t_symb), deriv, modules=['mpmath'])

def expand_bracket(func, cond, args, start, mult, maxiter):
    x = start
    if cond(func(x, *args)):
        raise Exception('invalid start point')
    for _ in range(maxiter):
        last_x = x
        x *= mult
        if cond(func(x, *args)):
            return last_x, x
    raise Exception(f'Unable to find starting area for {func} with {args}')

def safe_underflow_wrapper(func):
    def inner(x, *args):
        y_mpf = func(x, *args)
        y_f = float(y_mpf)
        if y_f == 0 and y_mpf != 0:
            y_f = float(mpmath.sign(y_mpf)) * math.ulp(0)
        return y_f
    return inner

def root_after(function, derivative, t, lower_bound):
    t = np.clip(t, 1e-30, np.inf)
    if derivative(lower_bound, t) > 0:
        cond = lambda x: x <= 0
    else:
        cond = lambda x: x > 0
    lo, hi = expand_bracket(derivative, cond, (t,), lower_bound, 5, 100)
    rng = hi - lo
    xtol = rng * 1e-3
    res = scipy.optimize.bisect(
        safe_underflow_wrapper(derivative), lo, hi, args=(t,), xtol=xtol
    )
    return res

ka_d0 = expr_derivative(kern_a_expr(), w0c, r0c, n=0)
ka_d1 = expr_derivative(kern_a_expr(), w0c, r0c, n=1)
kb_d0 = expr_derivative(kern_b_expr(), w0c, r0c, n=0)
kb_d1 = expr_derivative(kern_b_expr(), w0c, r0c, n=1)

last_root_a = None
last_root_b = None

def pick_points_a(t):
    r_star = root_after(ka_d0, ka_d1, t, lower_bound=1e-15)
    return [r_star * 0.01, r_star * 0.9, r_star, r_star * 1.3, r_star * 20]

def pick_points_b(t):
    r_star = root_after(kb_d0, kb_d1, t, lower_bound=1e-15)
    return [r_star * 0.01, r_star * 0.03, r_star, r_star * 1.4, r_star * 20]

def hopfield_fn_auto(t, w0, r0, dens, amp, bias):
    r_cap = 1.0
    pts_a = pick_points_a(t)
    pts_b = pick_points_b(t)
    a_val, _ = quad(kern_a, 0, r_cap, args=(w0, r0, t), points=pts_a, epsabs=0, epsrel=1.49e-8)
    b_val, _ = quad(kern_b, 0, r_cap, args=(w0, r0, t), points=pts_b, epsabs=0, epsrel=1.49e-8)
    return amp * ((4 * np.pi * dens * a_val) * np.exp(4 * np.pi * dens * b_val)) + bias

dens = 3e-19
amp = 1
bias = 0

teor_t0 = dens * np.pi * w0c * r0c**3
calc_t0 = hopfield_fn_auto(0, w0c, r0c, dens, amp, bias)
print('Theoretical value for t=0:', teor_t0)
print('Calculated value for t=0:', calc_t0)

T = np.linspace(0, 100000000, 1000)
Y = np.vectorize(hopfield_fn_auto)(T, w0c, r0c, dens, amp, bias)
plt.plot(T, Y, color='r', linestyle='-', linewidth=0.7)
plt.yscale('log')
plt.xscale('log')
plt.show()

Почему это важно

Адаптивная квадратура мощна, но не всеведуща. Если не направлять семплер к масштабно-чувствительным областям или если допуски не соответствуют вычисляемым величинам, возникают кажущиеся произвольными зависимости от параметров отсечки. Понимание механики points и настроек точности позволяет сохранить физику и убрать артефакты.

Практическое резюме

Интегрируя ядра вроде тех, что в модели Томаса—Хопфилда, заставьте интегратор смотреть туда, где есть структура, — передайте опорные точки через points. Для очень малых величин ужмите epsabs до нуля и опирайтесь на относительную точность. Если нужен настоящий несобственный интеграл от 0 до бесконечности, разделите его на конечную часть с points и бесконечный хвост без них. Использовать np.vectorize можно; обычный цикл на Python столь же корректен и иногда прозрачнее. И наконец, когда в схожих ситуациях считаете выражения вида exp(x) - 1 при малых x, подумайте о np.expm1(x), чтобы избежать катастрофической потери значимости.

С этими корректировками графики становятся стабильными, проверка при t=0 сходится с теорией, и можно переходить к подгонке модели с уверенностью, а не спорить с интегратором.