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

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

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

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

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

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

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

$$\int \approx h \cdot \left( \frac{f(x_0)}{2} + f(x_1) + f(x_2) + \dots + f(x_{n-1}) + \frac{f(x_n)}{2} \right), \quad h = \frac{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 обязано быть чётным):

$$\int \approx \frac{h}{3} \cdot \left( f_0 + 4 f_1 + 2 f_2 + 4 f_3 + \dots + 4 f_{n-1} + f_n \right)$$

Веса чередуются: концы — $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 возвращает не только значение, но и оценку погрешности — и сам выбирает число узлов адаптивно.

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

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

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

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

Итог

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