Численное интегрирование функций с сингулярностями: метод Tanh-Sinh
Когда подынтегральная функция уходит в бесконечность
Классические методы численного интегрирования (квадратуры Ньютона-Котеса, формулы Симпсона или Гаусса) блестяще справляются с гладкими функциями, которые заданы на конечном отрезке. Однако в физических задачах математикам постоянно приходится сталкиваться с интегралами, которые имеют математические особенности — сингулярности. Например, при расчетах потенциала заряженной нити или в задачах механики разрушения (напряжения на кончике трещины) подынтегральная функция может уходить в бесконечность на краях отрезка или иметь вертикальные асимптоты. Классический алгоритм, столкнувшись с бесконечностью (делением на ноль), выдаст ошибку Not-a-Number (NaN) или будет сходиться катастрофически медленно.
Существует несколько фундаментальных путей обхода этой проблемы. Первый, аналитический подход — это выделение и вычитание особенности. Мы искусственно добавляем и вычитаем из подынтегральной функции простую сингулярную функцию, интеграл от которой нам заранее известен аналитически. В результате у нас остается разность двух функций, которая в точке сингулярности имеет предел, равный константе. Эта гладкая часть (регулярная) легко интегрируется численно, а аналитическая часть просто прибавляется к итоговому ответу. Но этот подход требует ручной математической работы с каждым отдельным интегралом.
Магия нелинейных замен переменных
Второй, гораздо более мощный алгоритмический подход основан на нелинейных заменах координат. Идея заключается в том, чтобы так исказить ось интегрирования, чтобы «растянуть» проблемные участки (около сингулярности) и свести бесконечности к нулю. Одной из первых успешных замен такого рода стала подстановка Игараси-Мори-Такахаси (IMT).
Но настоящим прорывом и шедевром японской вычислительной школы стало создание в 1974 году квадратурной формулы гиперболического тангенса-синуса (Tanh-Sinh quadrature, или метод двойной экспоненциальной интеграции). Этот метод был разработан Хидетаси Такахаси и Масатаке Мори. Идея заключается во введении сложнейшей подстановки: x = tanh(pi/2 * sinh(t)). Эта функция обладает потрясающим топологическим свойством: она отображает бесконечный интервал (от минус до плюс бесконечности) строго в конечный отрезок [-1, 1].
Двойное экспоненциальное затухание
При выполнении этой замены дифференциал интегрирования (якобиан) порождает функцию, которая убывает к нулю с феноменальной скоростью — по закону двойной экспоненты: exp(-C * exp(|t|)). Эта скорость убывания настолько сокрушительна, что она способна мгновенно подавить (помножить на ноль) абсолютно любую математическую сингулярность исходной функции на краях отрезка!
После такой замены переменной мы получаем новый, гладкий и колоколообразный интеграл на бесконечной прямой, для вычисления которого применяется самое простое правило трапеций (да, базовая формула трапеций на бесконечном интервале для гладких убывающих функций работает точнее сложных методов Гаусса). Алгоритм Tanh-Sinh позволяет с феноменальной легкостью брать интегралы с любыми концевыми сингулярностями, логарифмическими разрывами или бесконечными производными, достигая невероятной точности в тысячи знаков после запятой в системах произвольной разрядности. Сегодня этот алгоритм встроен в ядра таких мощных математических систем, как MPMath в Python и Boost.Math в C++.