Формула Симпсона: парабола вместо отрезка

Урок про формулу Симпсона — рабочую лошадку численного интегрирования с четвёртым порядком точности.

Формула Симпсона приближает интеграл, заменяя функцию на каждой паре соседних полосок параболой через три точки; составная формула имеет 4-й порядок точности.

Идея: ловить кривизну

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

Веса узлов в формуле — 1, 4, 2, 4, 2, ..., 4, 1 (краевые с весом 1, средние полосок — 4, стыки — 2), всё умножается на h/3. Число полосок n должно быть чётным (узлы идут парами).

import math

def симпсон(f, a, b, n):
    if n % 2:           # n должно быть чётным
        n += 1
    h = (b - a) / n
    s = f(a) + f(b)
    for i in range(1, n):
        s += (4 if i % 2 else 2) * f(a + i * h)   # нечётные*4, чётные*2
    return s * h / 3

# Проверка 1: интеграл sin от 0 до pi = 2
print("∫sin(0..π):", round(симпсон(math.sin, 0.0, math.pi, 10), 8), " истина = 2")
# Проверка 2: интеграл x^4 от 0 до 1 = 0.2 (степень 4 — уже не точно)
print("∫x⁴(0..1):", round(симпсон(lambda x: x**4, 0.0, 1.0, 10), 8), " истина = 0.2")

# Порядок: ошибка падает в ~16 раз при удвоении n
точно = 2.0
print("\nсходимость на ∫sin(0..π):")
for n in [4, 8, 16, 32]:
    приближение = симпсон(math.sin, 0.0, math.pi, n)
    print(f"  n={n:2}: ошибка = {abs(приближение - точно):.2e}")

Вывод:

∫sin(0..π): 2.00010952  истина = 2
∫x⁴(0..1): 0.20001333  истина = 0.2

сходимость на ∫sin(0..π):
  n= 4: ошибка = 4.56e-03
  n= 8: ошибка = 2.69e-04
  n=16: ошибка = 1.66e-05
  n=32: ошибка = 1.03e-06

При удвоении n ошибка падает примерно в 16 раз (4.6e-3 → 2.7e-4 → 1.7e-5) — это четвёртый порядок в действии. Сравните с трапециями из прошлого урока: при n=16 Симпсон даёт ошибку 1.7e-5, трапеции — 6.5e-4, почти в 40 раз хуже при той же сетке.

Откуда веса 1-4-2-4-1

Формула Симпсона на двух полосках выводится интегрированием параболы Лагранжа через три узла: ∫ ≈ (h/3)(f₀ + 4f₁ + f₂). В составной формуле такие тройки накладываются: средний узел каждой пары входит с весом 4, а узлы-стыки принадлежат двум соседним тройкам — поэтому их веса складываются в 2. Краевые узлы участвуют лишь раз — вес 1. Отсюда характерный узор 1, 4, 2, 4, ..., 4, 1.

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

Погрешность составного Симпсона — −(b−a)·h⁴·f⁽⁴⁾(ξ)/180: четвёртая производная и h⁴ объясняют и высокий порядок, и точность на кубиках (у них f⁽⁴⁾=0). Симпсон — частный случай более общей идеи квадратур Ньютона–Котеса: брать полином всё выше через равноотстоящие узлы. Но выше Симпсона эта идея буксует: формулы высокого порядка на равномерных узлах получают отрицательные веса и теряют устойчивость (родственник феномена Рунге). Поэтому на практике берут не «Симпсон 10-го порядка», а либо адаптивный Симпсон (следующий урок), либо квадратуры Гаусса с неравномерными узлами.

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

Нечётное число полосок. Симпсон требует чётного n (узлы идут тройками с перекрытием); иначе формула неверна.
Перепутать веса 4 и 2. Вес 4 — у середин полосок (нечётные индексы), 2 — у стыков (чётные внутренние); ошибка в индексации портит ответ.
Гнаться за порядком через Ньютона–Котеса высоких степеней. Они неустойчивы; лучше адаптивный Симпсон или Гаусс.

Итоги

Симпсон заменяет функцию параболой на каждой паре полосок — 4-й порядок (удвоение n → ошибка ÷16).
Точна для полиномов до 3-й степени; веса узлов — 1, 4, 2, 4, ..., 4, 1, множитель h/3, n чётно.
Это частный случай Ньютона–Котеса; выше Симпсона равномерные формулы становятся неустойчивыми.
На практике — основа адаптивного интегрирования; при той же сетке точнее трапеций на порядки.