Квадратуры: численные методы приближения интегралов

Есть класс определённых интегралов, которые не поддаются вычислению «вручную». Не потому, что вы делаете что-то не так. У некоторых функций, например e^(−x²) или sin(x)/x, вообще нет первообразной в замкнутом виде. В таких случаях нужен другой подход: вычислить функцию в нескольких тщательно выбранных точках и объединить эти значения, чтобы оценить площадь под графиком. Это и есть квадратура.

В этой статье разбираем, что такое квадратура, как работает общий каркас и основные методы — от интуитивного правила середины до удивительно точного гауссова подхода. Мы также обсудим составные методы, источники ошибок, адаптивные стратегии и то, как выбрать метод под конкретную задачу.

Что такое квадратура?

Слово «квадратура» — старое. Оно происходит от классической геометрической задачи построить квадрат той же площади, что и данная фигура, то есть «найти площадь». Сегодня этот термин означает численное приближение определённого интеграла.

Ключевая идея — подмена: вместо точной первообразной вы берёте значения функции f(x) в наборе точек, называемых узлами, умножаете каждое значение на соответствующий вес и суммируете результаты:

Термины квадратура и численное интегрирование по сути означают одно и то же. «Квадратура» чаще встречается в более математических контекстах; «численное интегрирование» — повседневный термин в научных вычислениях.

Как работает численная квадратура

Любое квадратурное правило действует на интервале [a, b]. Внутри этого интервала оно выбирает набор узлов xᵢ и назначает каждому вес wᵢ. Приближение — это взвешенная сумма выше. Различие между методами — где располагаются узлы и как назначаются веса. Вот и всё. Правило трапеций, правило Симпсона и гауссова квадратура — это вариации одного и того же шаблона.

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

Распространённые методы квадратур

Ниже методы расположены от самых интуитивных к более изощрённым.

Наглядное сравнение правил средней точки, трапеций и Симпсона. Изображение автора.

Правило средней точки

Возьмите середину [a, b], вычислите там f и умножьте на ширину (b − a). Предполагается, что значение функции в центре представляет среднее по всему интервалу.

Звучит почти слишком просто — и для одного интервала так и есть. Но вот первый сюрприз: правило средней точки точно для констант и линейных функций, а его ошибка масштаба O(h³) на интервал, где h = b − a. На гладких функциях это на самом деле лучше, чем правило трапеций, что часто удивляет.

Правило трапеций

Используем оба конца вместо центра. Соедините их прямой и посчитайте площадь получившейся трапеции: (b − a) · [f(a) + f(b)] / 2. Точно для линейных функций, с ошибкой порядка O(h³) на интервал — того же порядка, что и у правила средней точки, хотя с большим коэффициентом.

Правило трапеций, вероятно, самый широко преподаваемый квадратурный метод, из-за чего его репутация несколько завышена. Для благоповодных функций оно неплохо. Но составное правило Симпсона почти всегда даёт лучшую точность при том же числе вычислений.

Правило Симпсона

Проведите квадратичный многочлен через три точки: два конца и середину. Проинтегрируйте этот многочлен точно:

Дополнительная оценка в центре стоит ещё одного вычисления функции, но отдача ощутима: правило Симпсона точно для многочленов до третьей степени, а ошибка на интервал масштаба O(h⁵) — на два порядка лучше, чем у правила трапеций. Для гладких функций это обычно первый метод, к которому стоит обратиться.

Гауссова квадратура

Здесь начинается самое интересное. Правила средней точки, трапеций и Симпсона используют равномерно расположенные узлы. Гауссова квадратура ставит другой вопрос: что если выбрать и расположение узлов, и веса так, чтобы максимизировать точность?

С n узлами гауссова квадратура (конкретно вариант Гаусса — Лежандра) способна точно интегрировать любой многочлен степени до 2n − 1. Два узла дают точность для кубических; пять узлов покрывают многочлены степени 9. Узлы — это корни многочленов Лежандра, и они сгущаются к концам интервала, а не располагаются равномерно. Дополнительная точность возникает благодаря свободному перемещению узлов, а не из-за того, где они «случайно» стоят. Их сгущение у концов — побочный эффект математики, и ниже мы увидим почему.

Вот как все четыре метода сравниваются:

Составные квадратурные методы

Одного применения правила на всём интервале [a, b] обычно недостаточно. Стандартное решение: разбить интервал на n равных подынтервалов шириной h = (b − a)/n и применить базовое правило на каждом из них.

Составные правила средней точки и трапеций

Для составного правила трапеций вычислите f в n+1 равномерно расположенных точках и просуммируйте n трапеций. Суммарная ошибка падает до O(h²). Уменьшите h вдвое — ошибка уменьшится вчетверо. Составное правило средней точки следует той же логике и имеет тот же порядок глобальной ошибки.

Составное правило Симпсона

Составной вариант требует чётного числа подынтервалов. Вычислите f в n+1 точках и получите глобальную ошибку порядка O(h⁴). Деление h пополам уменьшает ошибку в 16 раз. Поэтому составное правило Симпсона часто выбирают по умолчанию, когда нужна умеренная точность на гладкой функции.

Приведённый ниже код на Python показывает, как приближение улучшается при увеличении числа подынтервалов для ∫[0 до π] sin(x) dx (точное значение: 2.0):

import numpy as np

def composite_simpson(func, lower, upper, num_intervals):
    """Approximate integral using composite Simpson's rule."""
    if num_intervals % 2 != 0:
        raise ValueError("num_intervals must be even")
    nodes = np.linspace(lower, upper, num_intervals + 1)
    values = func(nodes)
    h = (upper - lower) / num_intervals
    weights = np.ones(num_intervals + 1)
    weights[1:-1:2] = 4   # odd-indexed interior points
    weights[2:-2:2] = 2   # even-indexed interior points
    return (h / 3) * np.dot(weights, values)

exact = 2.0
for n in [4, 8, 16, 32]:
    approx = composite_simpson(np.sin, 0, np.pi, n)
    print(f"n={n:3d}: approx={approx:.10f}, error={abs(approx - exact):.2e}")

Всего с 32 интервалами ошибка падает ниже 10⁻⁶. Это сходимость порядка O(h⁴).

Гауссова квадратура

Идея гауссовой квадратуры заслуживает большего, чем строчка в таблице.

Почему равномерные узлы не всегда оптимальны

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

Для n узлов правило Гаусса — Лежандра выбирает n корней n-го многочлена Лежандра Pₙ(x) как узлы, затем вычисляет веса так, чтобы любой многочлен степени ≤ 2n − 1 интегрировался без ошибки. На практике numpy.polynomial.legendre.leggauss(n) возвращает и то и другое. Узлы сгущаются к концам интервала [−1, 1], стандартного опорного интервала, где ошибка полиномиальной интерполяции обычно наибольшая и где повышенная плотность выборки даёт максимальную отдачу.

Когда использовать

Гауссова квадратура лучшая, когда у вас гладкая аналитическая функция и нужна высокая точность при малом числе вычислений. Пять точек Гаусса — Лежандра могут превзойти составное правило Симпсона с десятками интервалов на «хорошо ведущих себя» интеграндах. Но она плохо подходит для функций с изломами, разрывами или когда данные доступны только в фиксированных равномерно расположенных точках. В таких случаях тщательно выбранное расположение узлов играет против вас.

Ошибка и точность квадратуры

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

Источники ошибок включают:

• Дискретизационная ошибка (её также называют усечением) — это разрыв между точным интегралом и тем, что вернёт ваше квадратурное правило в идеальной арифметике. Именно её описывают формулы O(hᵏ). Она уменьшается по мере добавления точек или сужения интервалов.
• Ошибка плавающей запятой — это другая проблема. Суммирование множества малых величин в арифметике конечной точности вносит округления, которые не убывают с h и могут даже расти. В составных правилах с тысячами подынтервалов вы можете упереться в «пол», где шум округления доминирует и добавление интервалов перестаёт помогать. Для большинства научных задач вы достигнете нужной точности задолго до этого предела. Но он существует.

Гладкость функции, ширина интервала и выбор метода влияют на дискретизационную ошибку. «Грубая» функция с резкими особенностями требует больше интервалов или адаптивного подхода.

Адаптивная квадратура

Равномерная сетка тратит вычисления впустую, когда интегранд почти везде гладкий, но сложен в нескольких местах.

Как это работает

Примените квадратурное правило к интервалу, затем примените более высокого порядка (или то же правило на двух половинах) и сравните. Разница оценивает локальную ошибку. Если она ниже допуска, интервал принимается. Иначе разделите его и рекурсивно обработайте каждую половину. Повторяйте, пока каждый подынтервал не удовлетворит допуску.

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

scipy.integrate.quad()

В Python scipy.integrate.quad() реализует адаптивную гауссову квадратуру (конкретно алгоритм DQAGSE из QUADPACK). Передайте вашу функцию и пределы; получите оценку и границу ошибки. Для большинства одномерных интегралов это правильный инструмент — без ручной настройки.

from scipy import integrate
import numpy as np

# Integrate a function with a peak near x=0.5
def peaked_func(x):
    return np.exp(-100 * (x - 0.5)**2)

result, error_bound = integrate.quad(peaked_func, 0, 1)
print(f"Estimate: {result:.10f}")
print(f"Error bound: {error_bound:.2e}")

quad() уверенно справляется с этим узким пиком. Равномерному правилу Симпсона потребовались бы сотни подынтервалов, чтобы сопоставиться по точности.

Как выбрать квадратурный метод

Для одномерных интегралов в Python начинайте с scipy.integrate.quad(), если нет веских причин иначе. Это адаптивно, точно и хорошо проверено. Всё ниже актуально, когда вы реализуете что-то сами, работаете на языке без хорошей адаптивной библиотеки или оптимизируете под частный случай.

• Гладкие аналитические функции с известной формулой: гауссова квадратура даёт наибольшую точность на вычисление. Если нужна умеренная точность и простота, составное правило Симпсона — надёжный выбор по умолчанию.
• Функции с изломами или разрывами: методы с равномерными узлами (трапеции, Симпсон) часто превосходят гауссову квадратуру, поскольку фиксированное расположение узлов не «ломается» из-за нерегулярности. Разбейте интервал интегрирования в точках разрыва и обработайте части отдельно.
• Данные в фиксированных равномерных точках: гауссова квадратура не подходит. Вы не можете вычислять функцию в произвольных узлах, поэтому используйте составное правило трапеций или Симпсона.
• Интегралы повышенной размерности: квадратуры быстро упираются в «проклятие размерности». При n точках на измерение и d измерениях нужно nᵈ вычислений. Десять измерений по десять точек — это 10¹⁰ вычислений, что нереалистично. Монте-Карло масштабируется гораздо лучше выше примерно пяти измерений, там и происходит переключение.

Квадратура и родственные численные методы

Квадратура не существует в вакууме. Полезно видеть, где она находится относительно других инструментов, которые вы, вероятно, рассматриваете.

Квадратура vs аналитическое интегрирование

Аналитическое интегрирование находит точную первообразную; квадратура приближает число. Если точная первообразная существует и вы можете её найти, используйте её. Квадратура уместна, когда точное интегрирование невозможно или нецелесообразно.

Квадратура vs интегрирование Монте-Карло

Оба приближают определённые интегралы, но по-разному. Квадратура детерминирована: она вычисляет f в фиксированных заранее выбранных узлах, а сходимость связана с гладкостью функции и числом узлов. Монте-Карло берёт случайные выборки и усредняет их, причём ошибка убывает как O(1/√N) независимо от размерности. Это медленнее в малых размерностях, но не зависит от размерности, что и делает метод практичным примерно выше пяти измерений.

Квадратура vs интерполяция

Некоторые квадратурные правила, включая правило Симпсона, работают, подбирая многочлен по точкам вычислений и интегрируя его точно. Но цель отличается от интерполяции. Интерполяция восстанавливает функцию; квадратура нужна только ради числа. Многочлен — это инструмент, а не результат.

Заключение

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

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

В Python scipy.integrate.quad() — правильный выбор по умолчанию для одномерных задач. Всё остальное в этой статье — для понимания того, что происходит «под капотом», или когда вы реализуете что-то сами. Если хотите продолжить развивать численные навыки с SciPy, пройдите наш курс Introduction to Optimization in Python.