Математика / Численное интегрирование

Квадратура tanh–sinh: двойная экспоненциальная замена

Квадратура tanh–sinh переносит отрезок (−1,1) на всю прямую и применяет трапеции к преобразованной функции. Узлы быстро сгущаются к концам, помогая работать с интегрируемыми концевыми особенностями.

Опубликовано: Обновлено:

Формула

$$\int_{-1}^{1}f(x)\,dx\approx h\sum_{k=-N}^{N}f(\phi(kh))\phi^{\prime}(kh),\quad\phi(t)=\tanh\!\left(\frac\pi2\sinh t\right)$$
Схема Влияние длины суммы при h=0,5
00.250.2486N=10.0186N=22.60e−4N=36.50e−6N=4

Модули ошибок для f=1 на [−1,1]: N=1,2,3,4. Последние столбцы быстро уменьшаются.

Обозначения

$φ(t)$
монотонная замена из вещественной прямой в (−1,1)
$h$
шаг равномерной сетки в переменной t
$N$
число шагов в каждую сторону от нуля
$φ′(t)$
якобиан: (π/2)cosh(t)/cosh²((π/2)sinh(t))

Условия применения

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

Ограничения

Нельзя переносить утверждения о быстрой сходимости на произвольные разрывы внутри интервала. Численное округление φ(t) до±1 может вызвать бесконечное значение f ещё до умножения на малый якобиан. Нужно раздельно контролировать шаг h, ширину Nh и переполнение гиперболических функций.

Подробное объяснение

Замена φ(t)=tanh((π/2)sinh t) переводит t=0 в x=0, а бесконечные концы прямой — в±1. Производная всегда положительна. Поэтому исходный интеграл превращается в интеграл f(φ(t))φ′(t) по всей прямой, к которому применяют равномерную сумму с шагом h.

При больших |t| расстояние от φ(t) до края становится чрезвычайно малым, а якобиан быстро убывает. Именно сочетание двух гиперболических функций создаёт двойной экспоненциальный масштаб. Узлы на исходном отрезке оказываются плотнее возле концов, хотя сетка по t остаётся равномерной.

В сумме есть две ошибки. Ограничение индексами от−N доN отбрасывает хвост преобразованного интеграла. Конечный шаг h даёт ошибку квадратуры на оставшейся прямой. Если хвост уже мал, увеличение N при фиксированном грубом h почти ничего не исправит: требуется уменьшить шаг.

Это видно даже на f=1, чей интеграл равен 2. При h=1 сумма с N=1 даёт около 2,03084, а при N=2 —около 2,03137. Добавление точек не приблизило результат к 2: доминирует сеточная ошибка. При h=0,5 и N=4 получается 2,00000650, что значительно точнее.

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

Как пользоваться формулой

  1. Приведите исходный конечный отрезок к [−1,1] и выпишите, как изменятся функция и якобиан при линейной замене координаты.
  2. Выберите шаг h и индекс усечения N; для каждого k вычислите t=kh, затем φ(t) и точную формулу её производной.
  3. Сложите h·f(φ(t))·φ′(t) по симметричным индексам; не вводите половинные веса на искусственных концах усечённой суммы.
  4. Проверьте усечение увеличением N, затем отдельно уменьшите h; следите за округлением узлов до±1 и переполнением cosh.

Историческая справка

Хидэтоси Такахаси и Масатакэ Мори предложили двойную экспоненциальную замену для численного интегрирования в 1974 году. В обзоре Мори «Discovery of the Double Exponential Transformation and Its Developments», опубликованном в 2005 году, изложены происхождение подхода и его последующее развитие.

Основная идея соединяет удачную замену переменной с эффективностью равномерной формулы трапеций на всей прямой. Исторический акцент здесь относится к выбору преобразования, которое меняет распределение узлов и поведение хвостов. Современная запись tanh–sinh сохраняет эту структуру. Она не означает, что любое преобразование с гиперболическими функциями автоматически даёт ту же скорость сходимости или численную устойчивость.

Историческая линия формулы

Авторство двойной экспоненциальной формулы связывают с Такахаси и Мори и их работой 1974 года. Обзор Мори 2005 года является прямым свидетельством истории метода. Термин tanh–sinh описывает конкретную композицию функций, а не общее название всех квадратур с заменой переменной.

Пример

Проверим квадратуру на f(x)=1, I=2. Берём h=0,5, N=4, то есть девять значений k от−4 до 4. Значения f равны единице, поэтому остаётся сумма производных замены. Формула расчёта: I_h=0,5·Σ_{k=−4}^{4}[(π/2)cosh(0,5k)/cosh²((π/2)sinh(0,5k))]. Численно получаем 2,0000065048. Ответ — приближение 2, абсолютная ошибка около 6,50·10⁻⁶. Проверка малого индекса: при том же шаге и N=3 выходит 1,9997403043; добавление крайних членов уменьшает ошибку усечения. Однако при h=1 увеличение N с 1 до 2 не улучшает близость к 2. Значит, контроль одной только длины суммы недостаточен: шаг и усечение являются независимыми параметрами точности.

Частая ошибка

Не забывайте φ′: без якобиана вычисляется другая сумма, не приближающая исходный интеграл. Не увеличивайте только N, если ошибка определяется грубым шагом. Не вычисляйте особую функцию в округлённом узле x=1 с последующим умножением на нуль. Не переносите свойства концевых особенностей на произвольный разрыв внутри отрезка.

Практика

Задачи с решением

Центральный член

Условие. Найдите вклад k=0 при f=1 и h=0,5.

Решение. φ(0)=0, φ′(0)=π/2. Вклад равен hφ′(0)=π/4≈0,7853981634.

Ответ. π/4.

Нечётный множитель

Условие. Что даст любая симметричная усечённая сумма для f(x)=x?

Решение. φ нечётна, φ′ чётна. Вклады k и−k сокращаются, центральный член нулевой. Сумма равна 0, как и точный интеграл.

Ответ. 0 при любых h>0 и целых N≥0 в точной арифметике.

Дополнительные источники

  • https://www.kurims.kyoto-u.ac.jp/~prims/pdf/41-4/41-4-38.pdf
  • https://people.maths.ox.ac.uk/trefethen/sirev56-3_385.pdf

Связанные формулы

Математика

Составная формула трапеций для вычисления интеграла

$T_n=h\left[\frac{f(a)+f(b)}2+\sum_{i=1}^{n-1}f(a+ih)\right],\quad h=\frac{b-a}{n}$

Составная формула трапеций приближает определённый интеграл суммой площадей под линейными отрезками между соседними узлами. Концевые значения получают половинный вес.

Математика

Адаптивный Симпсон: локальная оценка ошибки

$S_1=S(a,b),\ S_2=S(a,m)+S(m,b),\quad \varepsilon_{\mathrm{est}}=\frac{|S_2-S_1|}{15},\quad I\approx S_2+\frac{S_2-S_1}{15}$

Адаптивный метод Симпсона сравнивает одну параболическую оценку с двумя оценками на половинах. Разность даёт локальный индикатор ошибки, по которому алгоритм уточняет только те участки, где не достигнут допуск.

Математика

Квадратура Кленшоу—Кёртиса по пяти косинусным узлам

$\int_{-1}^{1}f(x)\,dx\approx\frac{f(-1)+f(1)}{15}+\frac8{15}\left[f(-1/\sqrt2)+f(1/\sqrt2)\right]+\frac45f(0)$

Квадратура Кленшоу—Кёртиса интегрирует интерполяционный полином на косинусной сетке. Пятиузловая запись включает концы отрезка и точно работает для полиномов степени до пяти.