Численное интегрирование: трапеции и метод Симпсона

Площадь под кривой, когда первообразной нет под рукой: считаем интеграл численно — трапециями и по Симпсону.

Численное интегрирование (квадратура) — приближённое вычисление определённого интеграла суммированием площадей простых фигур под графиком.

Зачем считать интеграл численно

Аналитически взять интеграл получается далеко не всегда: у функции e^(-x²) (колокол Гаусса) вообще нет первообразной в элементарных функциях. А площадь под ней нужна постоянно — в статистике, физике, инженерии. Численные методы дают ответ для любой функции, которую можно вычислить в точках.

Метод трапеций

Разобьём отрезок [a, b] на n кусочков. На каждом заменим кривую прямой — получим трапецию, площадь которой легко считается. Сумма площадей трапеций ≈ интеграл:

∫ ≈ h * ( f(x0)/2 + f(x1) + f(x2) + ... + f(x_{n-1}) + f(xn)/2 ),   h = (b-a)/n

Концы берутся с весом 1/2, внутренние точки — с весом 1.

def f(x):
    return x * x          # интеграл x^2 от 0 до 1 = 1/3 = 0.3333...

a, b, n = 0.0, 1.0, 4
h = (b - a) / n
total = 0.5 * (f(a) + f(b))   # концы с весом 1/2
for i in range(1, n):
    total += f(a + i * h)     # внутренние с весом 1
print("Трапеции (n=4):", round(total * h, 6))
print("Точно         :", round(1/3, 6))

Вывод:

Трапеции (n=4): 0.34375
Точно         : 0.333333

Метод Симпсона: парабола вместо прямой

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

∫ ≈ (h/3) * ( f0 + 4*f1 + 2*f2 + 4*f3 + ... + 4*f_{n-1} + fn )

Веса чередуются: концы — 1, нечётные узлы — 4, чётные — 2.

def f(x):
    return x * x

a, b, n = 0.0, 1.0, 4      # n чётное
h = (b - a) / n
total = f(a) + f(b)
for i in range(1, n):
    coef = 4 if i % 2 == 1 else 2
    total += coef * f(a + i * h)
result = total * h / 3
print("Симпсон (n=4):", round(result, 6))
print("Точно        :", round(1/3, 6))

Вывод:

Симпсон (n=4): 0.333333
Точно        : 0.333333

Для многочлена 2-й степени Симпсон даёт точный ответ даже при n=4 — потому что параболой парабола приближается идеально.

А в SciPy — одна строка

from scipy.integrate import quad
import numpy as np

value, error = quad(lambda x: x**2, 0, 1)
print(value)   # 0.3333333333333333
print(error)   # ~3.7e-15 — оценка погрешности

# даже для функции без первообразной:
g, _ = quad(lambda x: np.exp(-x**2), 0, 1)
print(g)       # 0.7468241328124271

quad возвращает не только значение, но и оценку погрешности — и сам выбирает число узлов адаптивно.

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

Точность метода зависит от того, многочленом какой степени он «думает» функцию. Трапеции точны для линейных функций (ошибка ~ h²), Симпсон — для кубических (ошибка ~ h⁴). Удвоив число узлов, трапеции уменьшают ошибку в 4 раза, а Симпсон — в 16. Поэтому Симпсон при той же работе куда точнее. Промышленный quad идёт ещё дальше: использует адаптивную квадратуру Гаусса-Кронрода (QUADPACK), которая сама сгущает узлы там, где функция «капризничает», и оценивает собственную погрешность.

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

Нечётное n в Симпсоне. Метод требует чётное число интервалов — иначе формула неверна.
Особенности и разрывы. Если у функции внутри интервала разрыв или бесконечность, обычные квадратуры врут; quad умеет работать с особыми точками, но ему надо подсказать.
Слишком мало узлов. n=4 хорошо для гладких функций; для «дёрганых» нужно много больше.

Итог

Численное интегрирование суммирует площади простых фигур под графиком.
Трапеции (прямые, ошибка ~ h²) проще, Симпсон (параболы, ошибка ~ h⁴) точнее.
Симпсон требует чётного n и точен для многочленов до 3-й степени.
В SciPy — quad: адаптивный, с оценкой погрешности, работает даже без первообразной.