Численное интегрирование: прямоугольники, трапеции, Симпсон

Датчик расхода топлива пишет в бортовой журнал отсчёты каждую секунду: 3,2 л/ч, 3,1, 2,9… Сколько топлива сожжено за сутки? Это интеграл. Но функции нет — есть столбик из 86 400 чисел, и формула Ньютона–Лейбница беспомощна. Второй случай хуже: функция есть, а первообразной — нет. Колокол $e^{-x^2}$ из теории вероятностей, $\frac{\sin x}{x}$ из теории колебаний: их интегралы не выражаются через элементарные функции никак, ни за какое время стараний.

Почему формулы Ньютона–Лейбница недостаточно#

Признаемся честно: если $\int e^{-x^2}dx$ вам удалось взять, список элементарных функций придётся пополнить. Первообразной у гауссова колокола в этом классе не существует — это доказанный факт, а не недостаток фантазии. При этом значение ${\int }_0^1{e}^{-x^2}dx$ существует, и его умение вычислить нужно каждой программе, которая считает вероятности.

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

Прямоугольники: три варианта и три точности#

Отрезок $[a,b]$ делим на $n$ частей шага $h=\frac{b-a}{n}$ с узлами $x_i=a+ih$. Высоту каждого прямоугольника берём в одной из трёх точек — отсюда и три метода:

Левые и правые прямоугольники имеют точность $O(h)$: при уменьшении шага вдвое ошибка падает лишь вдвое, а главный член ошибки пропорционален $f^{\prime }(\xi )$. Средние прямоугольники точнее на целый порядок — $O(h^2)$. Причина красивая: разложение по формуле Тейлора в центре отрезка даёт $f(x_0+t)=f(x_0)+t{f}^{\prime }(x_0)+\dots$, а линейный член при интегрировании по симметричному отрезку $[-h/2,h/2]$ сокращается. Сдвиг точки выборки в середину бесплатно убивает главный источник ошибки — остаётся квадратичный.

Проверим всё это на конкретном интеграле ${\int }_1^2\frac{dx}{x}=\ln 2\approx 0,693147$ при $n=4$ ($h=0,25$). Все суммы здесь считаются вручную — пять значений функции $1/x$: $1$; $0,8$; $0,6667$; $0,5714$; $0,5$.

Тот же набор методов на ${\int }_0^1{e}^{x}dx=e-1\approx 1,718282$ при $n=4$ ($h=0,25$) — здесь видны и знаки ошибок, и их масштаб:

Знаки оправдали теорию: возрастающая и выпуклая экспонента занижается левыми и средними прямоугольниками и завышается правыми и трапециями. Средние точнее трапеций вдвое — отношение $8,9/4,5$ близко к двум. Симпсон отрывается на два с лишним порядка, хотя расходует те же пять значений функции, что и все остальные методы таблицы.

Метод трапеций: откуда берётся формула ошибки#

Трапеции соединяют точки $(x_i,f(x_i))$ и $(x_{i+1},f(x_{i+1}))$ хордой и складывают площади получившихся трапеций:

Формула погрешности выводится в две руки с рядом Тейлора. Растянем один отрезок:

1. 1
   Разложим $f(x_{i+1})$ около $x_i$: $f(x_i)+h{f}^{\prime }(x_i)+\frac{h^2}{2}{f}^{\prime \prime }(x_i)+\frac{h^3}{6}{f}^{\prime \prime \prime }(x_i)+\dots$
2. 2
   Точная площадь: ${\int }_{x_i}^{x_i+h}fdx=hf(x_i)+\frac{h^2}{2}{f}^{\prime }(x_i)+\frac{h^3}{6}{f}^{\prime \prime }(x_i)+\dots$ — это тот же ряд, проинтегрированный почленно
3. 3
   Площадь трапеции: $\frac{h}{2}(f(x_i)+f(x_{i+1}))=hf(x_i)+\frac{h^2}{2}{f}^{\prime }(x_i)+\frac{h^3}{4}{f}^{\prime \prime }(x_i)+\dots$
4. 4
   Разность: $\frac{h^3}{4}-\frac{h^3}{6}=\frac{h^3}{12}$ при $f^{\prime \prime }$ — трапеция завышает площадь на выпуклом участке
5. 5
   Складывая по всем отрезкам и сворачивая сумму вторых производных через интеграл, получаем итог: $T-I=\frac{(b-a)h^2}{12}{f}^{\prime \prime }(\xi )$

Итого трапеции — метод второго порядка: $O(h^2)$, то есть уменьшение шага вдвое снижает ошибку вчетверо. Числа это подтверждают: для ${\int }_1^{2}dx/x$ трапеции при $n=4$ ошиблись на $3,88\cdot 10^{-3}$, при $n=8$ — на $0,97\cdot 10^{-3}$. Отношение ровно четыре. Для экспоненты ${\int }_0^1{e}^{x}dx$: трапеции с $n=4$ дают $1,727222$ против точного $1,718282$ — ошибка $8,9\cdot 10^{-3}$.

Для измерений с неравномерным шагом трапеции всё равно работают: на каждом интервале площадь считается как $(x_{i+1}-x_i)\cdot \frac{f_i+f_{i+1}}{2}$ — половина суммы высот на длину интервала. Средним прямоугольникам нужны точки в серединах интервалов, которых в таблице может не оказаться; Симпсону — равномерная сетка и чётное число интервалов. Есть у трапеций и неожиданный талант: на периодических функциях, чей период укладывается в отрезок интегрирования целое число раз, они сходятся быстрее любого своего теоретического порядка — за это их любят в спектральных расчётах.

Формула Симпсона: парабола вместо хорды#

Хорда — многочлен первой степени. А что, если гнать через три соседние точки параболу? Симпсон так и делает: на паре отрезков $[x_0,x_2]$ строится интерполяционная парабола (тот же аппарат, что в интерполяции функций), интегрируется и даёт веса $1-4-1$ на трёх узлах. В составной формуле веса чередуются:

Число отрезков обязано быть чётным — иначе на последний отрезок параболы не хватит. Точность — четвёртый порядок, $O(h^4)$: шаг вдвое меньше — ошибка в 16 раз меньше. Бонус, который всегда удивляет: формула точна не только для парабол, но и для любых кубических многочленов — нечётные члены разложения при усреднении взаимно уничтожаются. Погрешность равна $-\frac{(b-a)h^4}{180}{f}^{(4)}(\xi )$.

Масштаб выигрыша виден на ${\int }_0^1{e}^{x}dx$ при $n=4$: Симпсон даёт $1,718319$ против точного $1,718282$ — ошибка $3,7\cdot 10^{-5}$. Трапеции на той же сетке ошиблись на $8,9\cdot 10^{-3}$: хуже в 240 раз при той же стоимости вычислений (пять значений функции у обоих).

Правило Рунге: как выбрать шаг автоматически#

Остался последний вопрос: сколько отрезков брать? Считать с запасом — жалко машинного времени, а в больших расчётах и стенки нет. Ответ дал Рунге: посчитайте интеграл дважды — с шагом $h$ и с шагом $h/2$ — и сравните. Для метода порядка $p$ ошибка более точного значения оценивается так:

Числа из примера с экспонентой: трапеции дали $T(0,25)=1,727222$ и $T(0,125)=1,720519$. Оценка: $\frac{|1,720519-1,727222|}{3}=0,00223$. Фактическая ошибка $T(0,125)$ равна $0,00224$ — оценка попала, не зная точного ответа. Заодно получен бонус — уточнение Ричардсона: $1,720519-0,00223=1,718285$ с ошибкой $3\cdot 10^{-6}$, лучше обоих исходных значений. Рекурсивное применение этой идеи — уже целая промышленная схема, а правила обращения с накопленной погрешностью разобраны в уроке про источники погрешностей.

Оговорка про честность: правило Рунге предполагает, что главный член ошибки действительно доминирует. Для разрывных, быстро осциллирующих или очень «крутых» функций разность $I_{h/2}-I_h$ может подводить, поэтому библиотеки после каждого удвоения $n$ сверяют новую оценку с предыдущей. Практический ориентир: если при последовательных удвоениях оценка стабильно падает примерно вчетверо (для трапеций) — вы в асимптотическом режиме, и цифрам можно верить. Оценка скачет — сетка ещё слишком грубая, продолжайте удвоение.

Какой метод выбрать под задачу#

Считайте не отрезки, а вычисления функции — именно они стоят денег. На ${\int }_0^1{e}^{x}dx$ Симпсон с восемью отрезками расходует девять значений функции и получает ошибку около $4\cdot 10^{-6}$; трапециям для той же точности понадобилось бы порядка двухсот значений. Когда каждое значение функции — минута счёта вспомогательной задачи, разница превращается из академической в суточную.

Порядок точности живёт на гладких функциях. Разрыв внутри отрезка роняет любой метод до первого порядка — режьте отрезок в точке разрыва и интегрируйте куски отдельно. Острые пики требуют сгущения сетки именно там, где они стоят; локальная проверка Рунге покажет, в каком месте ошибка не проходит. А если подынтегральная функция сама — решение дифференциального уравнения, интегрирование и решение ОДУ сливаются в одну задачу: об этом — финальный урок курса.

• Порядок точности важнее красивого названия: чтобы догнать Симпсона при $n=4$ на ${\int }_1^{2}dx/x$, трапециям нужно примерно $n=24$
• Ошибка имеет знак: выпуклая функция — трапеции завышают; знание знака позволяет строить вилки
• Разрыв или острый пик внутри отрезка обесценивает высокие порядки — делите отрезок на пике
• Для табличных данных доступна только та сетка, что измерена: трапеции — рабочая лошадка
• Погрешность считают по правилу Рунге — двумя расчётами, без обращения к точному ответу

Следующий логичный шаг — научиться по той же логике дифференцировать табличные функции: это численное дифференцирование с его ловушками шума. А если интеграл входит в дифференциальное уравнение, те же квадратуры превращаются в методы Рунге–Кутты — об этом наш урок про численное решение ОДУ. Потренировать ручные интегралы стоит в тренажёре интегралов, а термин интеграл в словаре пригодится, если определение подзабылось.