Адаптивное интегрирование и квадратуры Гаусса

Урок про два продвинутых подхода к интегрированию: адаптивное дробление по требованию и квадратуры Гаусса с оптимальными узлами.

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

Зачем адаптивность

Равномерная сетка расточительна: на гладких участках хватило бы пары узлов, а у резкого пика и сотни мало. Адаптивный метод решает это автоматически — он дробит отрезок только там, где нужно. Идея: посчитать интеграл на отрезке грубо и точнее (например, одним Симпсоном и Симпсоном на двух половинках). Если оценки близки — отрезок гладкий, принимаем результат. Если расходятся — функция сложна, рекурсивно делим отрезок пополам и применяем то же к половинкам. Так узлы сами сгущаются у пиков и разрежаются на ровном.

def адаптивный_симпсон(f, a, b, tol):
    def симпсон(a, b):
        c = (a + b) / 2
        return (b - a) / 6 * (f(a) + 4 * f(c) + f(b))
    def рек(a, b, целое, tol):
        c = (a + b) / 2
        левая = симпсон(a, c)
        правая = симпсон(c, b)
        if abs(левая + правая - целое) <= 15 * tol:   # критерий точности
            return левая + правая + (левая + правая - целое) / 15
        # иначе делим пополам и уточняем каждую половину
        return рек(a, c, левая, tol / 2) + рек(c, b, правая, tol / 2)
    return рек(a, b, симпсон(a, b), tol)

import math
# Интеграл 1/x от 1 до 2 = ln(2)
результат = адаптивный_симпсон(lambda x: 1.0 / x, 1.0, 2.0, 1e-10)
print(f"адаптивный:  {результат:.12f}")
print(f"ln(2)     :  {math.log(2):.12f}")
print(f"ошибка    :  {abs(результат - math.log(2)):.2e}")

Вывод:

адаптивный:  0.693147180560
ln(2)     :  0.693147180560
ошибка    :  8.33e-15

Адаптивный Симпсон достиг точности 8.3e-15 (почти машинный предел), потратив узлы экономно — больше там, где 1/x круче (у x=1), меньше у x=2. Слагаемое (левая+правая−целое)/15 — это экстраполяция Ричардсона (отдельный урок), бесплатно повышающая порядок.

Квадратуры Гаусса: умные узлы

Все формулы до сих пор брали узлы равномерно, подбирая только веса. Гаусс задал смелый вопрос: а что если оптимизировать и положения узлов? У формулы с k узлами тогда 2k свободных параметров (k весов + k позиций), и их хватает, чтобы точно интегрировать полиномы степени до 2k−1. Для сравнения: Симпсон с 3 узлами точен до степени 3, а Гаусс с 2 узлами — тоже до степени 3, но всего двумя вычислениями функции! Узлы Гаусса — это корни специальных ортогональных полиномов (Лежандра), они расположены неравномерно и сгущаются к краям.

import math

def гаусс2(f, a, b):
    # 2-узловая Гаусса-Лежандра: узлы ±1/√3 на [-1,1], веса = 1
    m = (b - a) / 2
    c = (a + b) / 2
    p = 1 / math.sqrt(3)
    return m * (f(c - m * p) + f(c + m * p))

# Всего ДВА вычисления функции, а кубику берёт точно:
print("∫x³(0..2):", round(гаусс2(lambda x: x**3, 0.0, 2.0), 8), " истина = 4")
print("∫eˣ(0..1):", round(гаусс2(math.exp, 0.0, 1.0), 8),
      " истина =", round(math.e - 1, 8))

Вывод:

∫x³(0..2): 4.0  истина = 4
∫eˣ(0..1): 1.71789638  истина = 1.71828183

Двумя вызовами функции Гаусс берёт x³ точно, а на eˣ ошибается лишь в 4-й цифре. С 5 узлами квадратура Гаусса даёт уже 9–10 верных цифр — недостижимая для трапеций экономия.

Как работает под капотом

Почему именно корни ортогональных полиномов? Если узлы — корни полинома Лежандра степени k, то ошибка квадратуры обнуляется для всех полиномов до степени 2k−1 благодаря свойству ортогональности. Веса при этом всегда положительны (в отличие от Ньютона–Котеса высоких степеней) — отсюда устойчивость. Минус Гаусса: при изменении k все узлы меняются, переиспользовать вычисления нельзя; это лечат вложенные правила (Гаусса–Кронрода), которые и стоят за scipy.integrate.quad — промышленным адаптивным интегратором.

Частые ошибки

• Равномерная сетка для функции с особенностью. Пик или почти-разрыв требует адаптивности; равномерная сетка либо неточна, либо чудовищно дорога.
• Думать, что «больше узлов = всегда лучше». Гаусс показывает: умно расставленные 2 узла бьют 10 равномерных; важно расположение, а не только число.
• Игнорировать критерий ошибки в адаптивном методе. Без сравнения грубой и точной оценок метод не знает, где дробить, и теряет смысл.

Итоги

• Адаптивное интегрирование дробит отрезок только там, где функция сложна, по критерию расхождения грубой и точной оценок.
• Квадратура Гаусса оптимизирует и узлы, и веса: k узлов точны до степени 2k−1.
• Узлы Гаусса — корни полиномов Лежандра, веса положительны (устойчивость).
• Промышленный интегратор (scipy.integrate.quad) сочетает адаптивность и Гаусса–Кронрода.