Квадратуры Ньютона

Содержание

Квадратуры Ньютона#

Задача численного интегрирования#

Требуется вычислить определенный интеграл

\[ I = \int_a^b f(x) dx \]

Интеграл предполагается собственным, то есть \(a, b \neq \infty\), \(|f(x)| < \infty\). Функция задана дискретно набором на сетке \(\{f(x_k) \}_{k=1}^N\).

Общий подход - правило квадратур#

Очевидно, что для численного приближения интеграла (вследствие линейности самого интеграла) общая формула для (численной) оценки будет

\[ Q^{(N)}=\sum_{k=1}^N w_k \cdot f\left(x_k\right), \]

Такая формула называется квадратурой. Таким образом для корректно-определенной квадратуры хотим стремление ошибки к нулю

\[ R^{(N)}=I-Q^{(N)} \rightarrow 0, \quad N \rightarrow \infty \]

Квадратуры Ньютона-Котеса#

Определим произвольную сетку на отрезке:

\[ a=x_0<x_1<\cdots<x_N=b \]

И разложим интеграл:

\[ I=\int_a^b f(x) d x=\sum_{k=1}^N I_k \]

\[ I_k=\int_{x_{k-1}}^{x_k} f(x) d x \]

Аппроксимируем каждый \(I_k\) как \(Q_k\) (теперь это интеграл по «маленькому» отрезку - гораздо более простая задача! Называется элементарная квадратура), и тогда итоговая составная квадратура будет

\[ Q^{(N)}=\sum_{k=1}^N Q_k \]

Осталось научиться аппроксимировать интегралы на «маленьком отрезке» - здесь нам поможет интерполяция.

Аппроксимация интегралов на малом отрезке. Элементарные квадратуры#

Если отрезок \([\alpha, \beta]\) достаточно мал, то можно приблизить нашу функции её интерполяционным многочленом (который интегрируется аналитически):

\[ \int_\alpha^\beta f(x) dx \approx \int_\alpha^\beta P(x) dx \]

Будем приближать функцию на отрезке \([\alpha,\beta]\) интерполяционным многочленом на равномерной сетке:

\(P_0(x) = f(\alpha)\). Получится формула прямоугольников: \(\displaystyle \int_\alpha^\beta f(x) dx \approx (\beta-\alpha) f(\alpha)\)
\(P_1(x) = f(\alpha) + (x-\alpha) \frac{f(\beta) - f(\alpha)}{\beta-\alpha}\). Получится формула трапеций: \(\displaystyle \int_\alpha^\beta f(x) dx \approx (\beta-\alpha) \frac{f(\alpha) + f(\beta)}{2}\)
Для интерполяционного многочлена второй степени (случай \(p = 2\)) получаем формулу Симпсона

\[\begin{split} \int_\alpha^\beta f(x) dx \approx (\beta-\alpha) \left[f(\alpha) w_0 + f\left(\frac{\alpha+\beta}{2}\right) w_1 + f(\beta) w_0\right],\\ w_0 = \frac{1}{2}\int_{-1}^1 \frac{x(x-1)}{2} dx = \frac{1}{6}\\ w_1 = \frac{1}{2}\int_{-1}^1 (1-x^2) dx = \frac{2}{3}\\ w_2 = \frac{1}{2}\int_{-1}^1 \frac{x(x+1)}{2} dx = \frac{1}{6} \end{split}\]

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

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

Пусть на отрезке \([a,b]\) задана равномерная сетка с интервалом \(h = \frac{b-a}{n}\). Применим к каждому отдельному интервалу формулу прямоугольников

\[ \int_a^b f(x) dx = \sum_{i=1}^n \int_{x_{i-1}}^{x_i} f(x) dx \approx h \sum_{i=1}^{n} f(x_{i-1}). \]

Применив к каждому формулу трапеций, получаем

\[\begin{split} \int_a^b f(x) dx = \sum_{i=1}^n \int_{x_{i-1}}^{x_i} f(x) dx \approx h \sum_{i=1}^{n} \frac{f(x_{i-1}) + f(x_i)}{2} =\\ = \frac{h}{2} \left(f(x_0) + 2f(x_1) + \dots + 2f(x_{n-1}) + f(x_n)\right) \end{split}\]

Объединив интервалы по два и, применив к каждой паре формулу Симпсона, получаем

\[\begin{split} \int_a^b f(x) dx = \sum_{i=1}^{n/2} \int_{x_{2i-2}}^{x_{2i}} f(x) dx \approx 2h \sum_{i=1}^{n/2} \frac{f(x_{2i-2}) + 4 f(x_{2i-1}) + f(x_{2i})}{6} =\\ = \frac{h}{3} \left(f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + \dots + 4f(x_{n-3}) + 2f(x_{n-2}) + 4f(x_{n-1}) + f(x_n)\right) \end{split}\]

Теоретическая погрешность составных формул#

Пример расчета погрешности#

Определим остаток следующей элементарной квадратуры:

\[\begin{split} \begin{aligned} R_k & =\int_{x_{k-1}}^{x_k} f(x) d x-h f\left(x_{k-1}\right) \\ & =\int_{x_{k-1}}^{x_k}\left[f(x)-f\left(x_{k-1}\right)\right] d x \end{aligned} \end{split}\]

Далее ряд Тейлора:

\[ f(x)=f\left(x_{k-1}\right)+f^{\prime}(\xi)\left(x-x_{k-1}\right), \quad \xi \in\left[x_{k-1}, x_k\right] \]

Итого,

\[\begin{split} \begin{aligned} R_k & =f^{\prime}(\xi) \int_{x_{k-1}}^{x_k}\left(x-x_{k-1}\right) d x \\ & =f^{\prime}(\xi) \int_0^h y d y=\frac{1}{2} f^{\prime}(\xi) h^2 \end{aligned} \end{split}\]

Пусть \(M_1=\max _{x \in[a, b]}\left|f^{\prime}(x)\right|\). Тогда погрешность составной квадратуры будет

\[\begin{split} \begin{aligned} \left|R^{(N)}\right|=\left|\sum_{k=1}^N R_k\right| & \leqslant \frac{1}{2} M_1 h^2 N \\ & =\frac{1}{2} M_1(b-a) h \end{aligned} \end{split}\]

Несложно проделать эту же процедуру и для других методов:

\(p\)	\(n\)	Название	Погрешность
\(0\)	\(\forall\)	Прямоугольников	\(\frac{(b-a) h}{2} M_1\)
\(1\)	\(\forall\)	Трапеций	\(\frac{(b-a) h^2}{12} M_2\)
\(2\)	\(n = 2m\)	Симпсона	\(\frac{(b-a) h^4}{180} M_4\)
\(3\)	\(n = 3m\)	Формула 3/8	\(\frac{(b-a) h^4}{80} M_4\)

Показатель степени \(h\), как обычно, называется порядком сходимости метода.

Точность квадратурных формул#

Алгебраической степенью точности квадратурной формулы называют такое число \(d\), что квадратурная формула точна для всех многочленов степени \(d\), но для некоторых многочленов степени \(d+1\) уже не точна. Из вида остаточного члена заключаем, что метод прямоугольников имеет \(d=0\), метод трапеций \(d=1\), методы Симпсона и 3/8 \(d=3\).

Графики фактических погрешностей для простых квадратур#

def rectangle(f, h):
    return h * sum(f[:-1])

def trapezoid(f, h):
    return 0.5 * h * (f[0] + 2 * sum(f[1:-1]) + f[-1])

def simpson(f, h):
    assert (len(f)-1) % 2 == 0
    return h/3. * (f[0] + 4 * sum(f[1:-1:2]) + \
                   2 * sum(f[2:-2:2]) + f[-1])

def threeeights(f, h):
    assert (len(f)-1) % 3 == 0
    return 3*h/8. * (f[0] + 3 * sum(f[1:-1:3]) + \
                3 * sum(f[2:-1:3]) + 2 * sum(f[3:-3:3]) + f[-1])

../../_images/524d06507d2c7142dd1b748a4b44a1c82ba2c462a528c1328944c5b01d2f2926.png