Численное интегрирование

Землю нельзя положить на весы. Но сейсмологи по тому, как сквозь неё проходят волны землетрясений, знают, какая у неё плотность на каждой глубине: в 1981 году Адам Дзевоньский и Дон Андерсон свели эти измерения в модель PREM — десять слоёв, от железного ядра до океана. Масса — это сумма тонких сферических оболочек: интеграл 4πr²ρ(r) от центра до поверхности. Плотность в модели задана многочленами, поэтому точный ответ мы знаем и можем честно измерить ошибку каждого способа считать интеграл суммой. Выберите способ, двигайте число вычислений и переключите «каждый слой отдельно».

Начните с Симпсона и узлов, расставленных ровно от центра до поверхности. На 4097 вычислениях плотности он ошибается на 1,2·10⁻⁴ — почти как прямоугольники, у которых 1,1·10⁻⁴, — и с ростом числа узлов ошибка падает медленно и неровно. Теперь включите «каждый слой отдельно»: та же формула на 1024 кусках ошибается на 5,6·10⁻¹³, а три точки Гаусса в каждом слое — 30 вычислений — дают все шестнадцать знаков. Формулы те же, данные те же; разница в том, знает ли метод, где у функции скачки. Об этом вся статья.

Идея

Отрезок режут на куски длиной h и на каждом заменяют функцию тем, что интегрируется точно: константой, отрезком прямой, параболой. Площади кусков складывают. Если функция гладкая, на маленьком куске она почти совпадает с заменой, и ошибка куска пропорциональна высокой степени h. Кусков 1/h, поэтому общая ошибка ведёт себя как h^p. Число p — порядок метода: чем он больше, тем быстрее точность растёт с числом вычислений.

Почему не взять первообразную? Потому что её часто нет. У гауссовой кривой e^−x²/2, у 1/ln x, у длины эллипса первообразных среди элементарных функций нет, а у функции, заданной таблицей замеров, нет даже формулы. Геометрия интеграла, формулы Ньютона — Лейбница и Симпсона подробно разобраны в главе об интеграле «Царицы наук»; здесь — что из этого выходит как из алгоритма: сколько стоит точность, как узнать ошибку, не зная ответа, и где формулы врут.

Три формулы и их порядок

Пусть отрезок [a, b] разбит на n кусков длиной h = (b − a)/n, и f_k — значение функции в точке a + kh.

• Прямоугольники (средних точек): на каждом куске высота — значение в его середине. Сумма h·(f(a + h/2) + f(a + 3h/2) + …).
• Трапеции: соседние точки соединяют отрезками. Сумма h·(f₀/2 + f₁ + … + f_n−1 + f_n/2).
• Симпсон: через каждые три соседние точки проводят параболу; n чётное, веса 1, 4, 2, 4, …, 2, 4, 1, умноженные на h/3.

Возьмём интеграл, которого нет в элементарных функциях: долю значений нормального распределения в пределах одного σ, ∫ от −1 до 1 плотности e^−x²/2/√(2π). Ответ 0,6826894921370859 известен по таблицам функции ошибок. Программа считает его тремя формулами и печатает ошибки:

   n    rectangles   trapezoids      Simpson   Runge for Simpson
   4    5.116e-03   -1.017e-02    3.686e-04     6.786e-04
   8    1.265e-03   -2.526e-03    2.148e-05     2.314e-05
  16    3.154e-04   -6.305e-04    1.320e-06     1.344e-06
  32    7.878e-05   -1.576e-04    8.216e-08     8.253e-08
  64    1.969e-05   -3.938e-05    5.130e-09     5.136e-09

Удвоили число кусков — ошибка прямоугольников и трапеций упала вчетверо, а Симпсона — в шестнадцать раз. Это и есть порядки 2, 2 и 4. Заметьте ещё одно: трапеции ошибаются вдвое сильнее прямоугольников и в другую сторону. Поэтому сумма двух третей прямоугольников и одной трети трапеций гасит главную ошибку, и получается ровно формула Симпсона с вдвое меньшим шагом. Откуда берётся h², видно из точной формулы ошибки.

Теорема говорит больше, чем «ошибка порядка h²». Она точно называет, чем платят: второй производной. Функции, у которых f″ ограничена, трапеции считают с ошибкой ~ h², а там, где производная взрывается или функция прыгает, гарантия пропадает. Для прямоугольников так же выводится ошибка +(b − a)h²f″/24 — вдвое меньше и с другим знаком, что мы и видели в таблице. Для Симпсона рассуждение длиннее (интегрировать по частям приходится четыре раза), а ответ −(b − a)h⁴f⁗(η)/180: платят четвёртой производной, и формула точна для любого кубического многочлена.

На графике внизу первой сцены по осям отложены логарифмы, и степенной закон «ошибка ≈ C/n^p» превращается в прямую с наклоном −p. Когда слои считаются отдельно, линии прямоугольников и трапеций идут вниз с наклоном −2, Симпсона — с наклоном −4, а Гаусс сразу лежит на дне, на уровне ошибок округления. Когда узлы расставлены вслепую, все четыре линии сбиваются в кучу с наклоном около −1. Причина — скачки плотности на границах слоёв: кусок, который перекрывает скачок высотой J, ошибается на величину порядка J·h, какую бы формулу на нём ни применили, а такой кусок на каждой границе есть всегда. Девять границ — девять ошибок порядка h, и порядок любой формулы падает до первого.

Ошибка без ответа: правило Рунге

В таблице выше ошибку легко посчитать: ответ известен. В настоящей задаче его нет, иначе незачем было бы считать. Выручает то, что ошибка метода устроена предсказуемо: если порядок p, то при шаге вдвое меньше она уменьшается в 2^p раз. Значит, два ответа с шагами h и h/2 отличаются в основном на ошибку, и по их разнице её можно восстановить.

Это правило Рунге: посчитай дважды, раздели разницу на 2^p − 1 — для трапеций на 3, для Симпсона на 15. Последний столбец таблицы выше — именно такая оценка ошибки Симпсона, без знания точного ответа: на 4 кусках она промахнулась почти вдвое, а начиная с 32 кусков совпадает с настоящей ошибкой с точностью до процента. Оговорка в условии важна: правило работает, когда h уже настолько мал, что ошибка ведёт себя как Ch^p. На скачке или на шумных данных этого режима нет, и оценка становится грубой прикидкой.

Оценку можно не только печатать, но и прибавить к ответу: I_h/2 + (I_h/2 − I_h)/3 для трапеций — это ровно формула Симпсона, и ошибка становится порядка h⁴. Повторяя приём с поправками всё более высоких порядков, получают метод Ромберга — таблицу, где каждый следующий столбец точнее предыдущего. Сам приём, уточнение по двум шагам, в общем виде описал Льюис Ричардсон в 1911 году; в учебниках его называют экстраполяцией Ричардсона.

Адаптивный шаг: дробить там, где трудно

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

Пройдите шаги на √x. Её производная у нуля бесконечна, и метод уходит к нулю всё глубже: полоса точек под осью густеет слева, а справа хватает нескольких кусков. На допуске 10⁻⁷ ему нужно 185 вычислений и ошибка 8·10⁻¹⁰; Симпсон с равным шагом на тех же 185 точках ошибается на 3·10⁻⁵. Переключитесь на Землю: метод ничего не знает о слоях, но сам находит все девять границ — около каждой он делит отрезок до самого дна, как бисекция, которая ищет корень.

Число 15 в условии остановки — из правила Рунге: у Симпсона p = 4, и ошибка двух половинок примерно в 2⁴ − 1 = 15 раз меньше их расхождения с целым. Разницу, делённую на 15, метод тут же прибавляет к ответу, как Ричардсон. Допуск при делении уменьшается вдвое, чтобы сумма ошибок всех кусков осталась в пределах исходного. На Земле вслепую он делает 1485 вычислений при допуске 10⁻⁶ и получает ошибку 1,3·10⁻⁹ — почти в сто тысяч раз меньше, чем Симпсон с равным шагом на 4097 точках. Глубина деления ограничена пятьюдесятью: на скачке расхождение не исчезает никогда, и без ограничения метод дробил бы кусок до размеров ошибки округления.

Первым рекурсивный адаптивный Симпсон напечатал Г. Ф. Кунцир в 1962 году — алгоритм 103 в журнале Communications of the ACM. Нынешний вид — с делением пополам и поправкой на одну пятнадцатую — придал ему Джеймс Лайнесс в 1969 году.

Ловушка: пик между узлами

Переключите сцену на узкий пик: парабола 1 + x² и на ней гауссов всплеск шириной 0,002 в точке 0,37. Метод смотрит в точки 0, 1/4, 1/2, 3/4 и 1, видит параболу, на которой формула Симпсона точна, — половинки сходятся с целым до последнего знака — и останавливается после пяти вычислений с ответом 1,3333333. Правильный ответ 1,3368782. Допуск хоть 10⁻¹²: ответ тот же. Подвигайте пик к 1/4 или 1/2 или сделайте его шире — тогда метод его заметит и раздробит.

Поэтому никакая квадратура не гарантирует точность без предположений о функции. Гарантия появляется, только если о функции что-то известно заранее: ограничена производная, нет пиков уже некоторой ширины, известны места скачков. Библиотеки так прямо и спрашивают: у scipy.integrate.quad есть параметр points, куда передают трудные места внутри отрезка — особенности и разрывы. Это тот самый переключатель «каждый слой отдельно» из первой сцены.

Узлы Гаусса: вдвое больше точности за те же вычисления

У всех формул выше узлы стоят через равные промежутки, и свобода есть только в весах. Формула по n точкам с удачными весами точна для многочленов степени до n − 1: веса подбирают так, чтобы интеграл совпал на 1, x, …, xⁿ⁻¹. Гаусс спросил: а если двигать и сами узлы? Тогда свободных чисел вдвое больше — n узлов и n весов, — и можно надеяться на точность для степени 2n − 1. Надежда сбывается, если поставить узлы в корни многочлена Лежандра.

Корни P_n программа находит методом Ньютона, стартуя из косинусов равномерной сетки углов, а сами P_n считает по трёхчленной формуле (k + 1)P_k+1 = (2k + 1)xP_k − kP_k−1. Три узла — 0 и ±√(3/5) с весами 8/9 и 5/9 — интегрируют x⁰, …, x⁵ точно и впервые ошибаются на x⁶. PREM в каждом слое — многочлен не выше третьей степени, после умножения на r² — не выше пятой. Поэтому три точки Гаусса на слой дают массу Земли точно, и десять слоёв — это тридцать вычислений.

Пример в конце программы взят из мемуара самого Гаусса 1815 года «Новый метод приближённого вычисления интегралов». Интеграл dx/ln x от 100 000 до 200 000 Гаусса интересовал не случайно: он примерно равен числу простых чисел в этом промежутке. Гаусс посчитал его по одной, двум, …, семи точкам вручную, и его ответ по семи точкам — 8406,2431211 — верен в девяти знаках. У нашей программы семь узлов дают ошибку 2,6·10⁻⁹, а Симпсон на тех же семи точках ошибается на 7·10⁻³. Простых чисел между 100 000 и 200 000 на самом деле 8392 — это уже вопрос о простых числах, а не о квадратурах.

nodes   Gauss            error        Simpson on as many nodes
    3   8406.2367752    -6.346e-03      4.421e-01
    5   8406.2431171    -3.779e-06      3.399e-02
    7   8406.2431208    -2.550e-09      7.061e-03

У формулы Гаусса есть неудобство: узлы для n и для 2n точек не совпадают, и при уточнении все вычисления приходится выбрасывать. В 1964 году Александр Кронрод нашёл, как дополнить n узлов Гаусса ещё n + 1 узлами, чтобы новая формула была точна до степени 3n + 1, а старые вычисления пошли в дело. Разница двух ответов — оценка ошибки почти бесплатно. На этой паре держится библиотека QUADPACK 1983 года, а через неё — scipy.integrate.quad: адаптивное деление, на каждом куске формула Кронрода по 21 точке.

Когда есть только таблица замеров

Всё, что было выше, требовало права вычислить функцию где угодно. С замерами так нельзя: точки уже стоят, где стояли, и ни адаптивный шаг, ни узлы Гаусса недоступны. Остаются трапеции и Симпсон, и вопрос меняется: не «какой формулой», а «хватает ли замеров».

Пример — отопление. Сколько тепла дом теряет за зиму, примерно пропорционально градусо-суткам: интегралу от разницы между 18 °C и уличной температурой по всем часам отопительного сезона. Ниже — Алматы с 15 октября 2024 по 15 апреля 2025, по часам. По всем 4368 часам выходит 3036 градусо-суток. А теперь представьте наблюдателя, который записывает температуру раз в сутки и считает, что она держалась так до следующей записи.

Запись в 14:00 недосчитывается почти 27 % тепла: днём теплее, чем в среднем за сутки, и наблюдатель каждый день видит самый тёплый час. Запись в 6 утра переоценивает на 19 %. А в 9 утра или в 7 вечера ошибка меньше процента, но это везение: в эти часы температура проходит через своё среднее за сутки. Посмотрите на строку под графиком: трапеции и Симпсон по тем же записям дают почти ту же ошибку. Формула высокого порядка не помогает, потому что ошибка здесь не от того, как соединены точки, а от того, когда их взяли: шаг равен периоду суточного хода температуры, и все записи попадают в одну и ту же его фазу.

Четыре записи в сутки, в 5, 11, 17 и 23 часа (это 0, 6, 12 и 18 часов всемирного времени — главные сроки наблюдений метеостанций), ошибаются на 0,7 %, восемь — на 0,2 %. Это неслучайно мало. Суточный ход повторяется каждые 24 часа, а сумма по N равноотстоящим точкам за период точна для постоянной и для любой гармоники, номер которой не делится на N: сумма e^2πikj/N по j от 0 до N − 1 равна нулю, как и интеграл. Четыре записи не ошибаются на волнах с периодом 24, 12 и 8 часов, а у суточного хода почти вся амплитуда в первых гармониках. Для гладкой периодической функции такая сумма сходится не как h², а быстрее любой степени; Ллойд Трефетен и Андре Вейдеман посвятили этому в 2014 году обзор «Экспоненциально сходящаяся формула трапеций». Симпсон на периодической функции даже хуже трапеций: его веса 4 и 2 портят равноправие точек.

Правило Рунге работает и на таблице: посчитайте интеграл по всем замерам и по каждому второму. Здесь это 3036,22 и 3036,24 градусо-суток — значит, часовых замеров с запасом хватает. Если же два ответа расходятся на проценты, никакая формула не спасёт: нужны замеры чаще. А если замеры идут неравномерно — датчик пропускал часы, — трапеции работают на любом шаге: каждая трапеция берёт свою ширину. Так устроены numpy.trapezoid и scipy.integrate.simpson, которым вместе со значениями передают и моменты замеров.

Код: Земля по слоям

Вот программа, которая стоит за первой сценой. Плотность — многочлены PREM, точный ответ — по первообразным каждого слоя, дальше четыре формулы вслепую и по слоям. Те же три точки Гаусса дают и всё остальное, что сейсмологи знают о недрах из одного профиля плотности.

moment of inertia I/(MR^2) = 0.3308
gravity: 9.82 m/s^2 at the surface, the most 10.69 at 3480 km from the centre
pressure in the centre: 364.1 GPa

Момент инерции 0,3308MR² — меньше, чем 0,4 у однородного шара: масса собрана к центру. Это число независимо измеряют по прецессии земной оси, и справочник NASA даёт тот же 0,3308. Сила тяжести в глубине считается по массе внутри радиуса r — снова интеграл, — и она удивляет: до границы ядра на глубине 2891 км она не падает, а растёт до 10,69 м/с². Оболочка над головой не тянет ни вверх, ни вниз, а масса под ногами убывает медленнее, чем растёт множитель 1/r²: почти треть массы Земли — в тяжёлом ядре. Давление в центре — вес всех слоёв над ним, интеграл ρ(r)g(r), где g сам интеграл: 364 гигапаскаля, три с половиной миллиона атмосфер. Интеграл внутри интеграла с теми же тремя точками на кусок: внутри слоя всё гладко, и 64 кусков на слой хватает с огромным запасом.

Сколько это стоит

Время любой квадратуры — это число вычислений функции: суммирование почти ничего не стоит, а функция бывает тяжёлой, вплоть до целой симуляции на каждую точку. Память — O(1) для формул с постоянным шагом и стек глубины деления для адаптивной. Главный вопрос поэтому: сколько вычислений нужно для точности ε.

На скачках и изломах все столбцы, кроме последнего, теряют смысл: порядок падает до первого, если не разрезать отрезок по особым точкам. Если вы знаете, где они, — режьте; не знаете — адаптивный метод найдёт их сам, хоть и ценой лишних вычислений.

В нескольких измерениях цена растёт быстро. Сетка по n точек на каждую из d осей — это n^d вычислений, и ошибка Симпсона в пересчёте на число точек N падает как N^−4/d. У метода Монте-Карло она падает как N^−1/2 в любой размерности, и при d = 8 эти два закона сравниваются. Дальше случайные точки выгоднее любой сетки.

Где это работает

• Научные библиотеки. scipy.integrate.quad — обёртка над QUADPACK: адаптивное деление, формула Гаусса — Кронрода на кусках, допуски по умолчанию около 1,5·10⁻⁸ и параметр points для разрывов.
• Таблицы замеров. numpy.trapezoid и scipy.integrate.simpson принимают моменты замеров вместе со значениями: энергия по показаниям мощности, путь по скорости, объём по расходу.
• Метод конечных элементов. Прочность моста или крыла считают, разбивая конструкцию на элементы; интегралы по каждому элементу берут не формулой, а в нескольких точках Гаусса.
• Фармакология. Сколько лекарства попало в кровь, оценивают площадью под кривой концентрации (AUC) по анализам, взятым через неравные промежутки, — трапециями.
• Строительство. Градусо-сутки отопительного периода задают, насколько тёплыми должны быть стены, — та же сумма, что в последней сцене, только по многолетним данным.

История

Трапеции старше интеграла почти на две тысячи лет. На глиняных клинописных табличках 350–50 годов до нашей эры астрономы вычисляли, насколько Юпитер сдвинется по небу за 60 дней, как площадь трапеции под графиком его суточной скорости, а момент, когда пройдена половина пути, находили, деля трапецию на две равные по площади. Расшифровал эти расчёты историк Матье Оссендрейвер в 2016 году. Параболу вместо прямой предложил Кеплер в 1615 году для объёма винных бочек; общую формулу напечатал Томас Симпсон в 1743-м — подробнее об этом в главе об интеграле «Царицы наук». Формулы с равными шагами до одиннадцати точек вывел Роджер Котс; его бумаги издали в 1722 году, после смерти автора, и с тех пор эти формулы называют формулами Ньютона — Котса.

Мемуар Гаусса о свободных узлах вышел в 1815 году, и вывел он их через цепные дроби. Вывод через ортогональные многочлены, по которому их учат теперь, дал Карл Якоби в 1826 году. Правило оценки ошибки по двум шагам предложил Карл Рунге в начале XX века, Льюис Ричардсон в 1911 году превратил его в способ уточнять ответ, а Вернер Ромберг в 1955 году — в таблицу повторных уточнений. Шестидесятые принесли компьютерные квадратуры: адаптивный Симпсон Кунцира в 1962 году, формулы Кронрода в 1964-м, а в 1983 году Роберт Писсенс, Элиза де Донкер-Капенга, Кристоф Юбергубер и Дэвид Каханер собрали библиотеку QUADPACK, которая работает внутри SciPy до сих пор.

И одна история о том, как медленно расходится знание. В 1994 году в журнале Diabetes Care вышла статья исследовательницы питания Мэри Тай: способ считать площадь под кривой сахара в крови, разбивая её на трапеции и складывая. Это была формула трапеций, о чём журнал вскоре напечатал письма с говорящими названиями вроде «Формула Тай — это формула трапеций». Тай ответила, что вывела способ сама в 1981 году, а напечатала потому, что коллеги уже пользовались им и хотели на что-то ссылаться. Ссылаются до сих пор: к 2025 году у статьи больше пятисот цитирований.

Задачи

В первой задаче настоящие часы Алматы с пропусками: интегрировать нужно по неравным шагам и аккуратно там, где температура пересекает 18 °C. Во второй — планета из слоёв: масса, момент инерции и давление в центре, где один интеграл стоит внутри другого.