Тема 17

Интерполационни квадратурни формули. Квадратурни формули на правоъгълниците, трапеците и Симпсон, представяне и оценка на грешките им.


страницата се нуждае от дописване/преглеждане


Бележки:
Формула на Лайбниц-Нютон:
$\displaystyle \int\limits_{a}^{b} f(x) dx = F(b) - F(a)$
където $F'(x) = f(x)$, т.е $F(x)$ е примитивна функция на $f(x)$.


В тази тема се разглежда въпросът за численото интегриране и се поставя следната интерполационна

Задача

Дадени са точките $a \leq x_1 < x_2 < \dots < x_n \leq b$ и съответните функционални стойности $\{f(x_i)\}_{i=1}^{n}$.
Търсим приближение за $I(f) = \int\limits_{a}^{b} f(x) dx$.

Решение:
Можем да построим $L_{n-1}(f;x) = \sum_{k=1}^n l_k(x) f(x_k)$ - интерполационният полином за $f(x)$ с възли $x_1, x_2, \dots, x_n$
Тук, припомням, $\displaystyle l_k(x) = \prod_{i = 1, i \neq k}^n \frac{x - x_i}{x_k - x_i}$

Тогава

(1)
\begin{align} \int\limits_{a}^{b} f(x) dx \approx \int\limits_{a}^{b} L_{n-1}(f;x) dx = \sum_{k=1}^{n} \int\limits_{a}^{b} l_k(x) dx . f(x_k) \end{align}

Получаваме

(2)
\begin{align} I(f) = \int\limits_{a}^{b} f(x) dx \approx \sum_{k=1}^n a_k f(x_k) \end{align}

където означаваме $\displaystyle \int\limits_{a}^{b} l_k(x) dx = a_k, \ k = 1, 2, \dots, n$

Дефиниция (квадратурна формула)

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

Формулата $Q(f) = \sum_{k=1}^n a_k f(x_k)$ е квадратурна формула

Точките $\{x_k\}_{k=1}^n$ се наричат възли на квадратурната формула.
Коефициентите $\{a_k\}_{k=1}^n$ се наричат коефициенти на квадратурната формула.

Дефиниция (интерполационна формула)

Ако квадратурната формула е $Q(f) = \sum_{k=1}^n a_k f(x_k)$, то тя се нарича интерполационна ако е получена посредством
$\displaystyle Q(f) = \int\limits_{a}^{b} L_{n-1}(f;x) dx$
т.е коефициентите и се получават по формулата
$\displaystyle a_k = \int\limits_{a}^{b} l_k(x) dx, \ k = 1, 2, \dots, n$

Дефиниция (точност на квадратурна формула)

Ще казваме, че квадратурната формула $Q$ е точна за $f$, ако остатъкът е равен на нула.

Остатъкът бележим с $R(Q;f) = I(f) - Q(f)$.

Теорема

Квадратурната формула $Q(f) = \sum_{k=1}^n a_k f(x_k)$ е интерполационна, т.с.т.к. тя е точна за $\forall f \in \Pi_{n-1}$

Представяне и оценка на грешката за интерполационни формули

От формулата на Нютон, която гласи следното
$f(x) = L_{n-1}(f;x) + f[x_1, x_2, \dots, x_n, x]. \omega(x)$,
чрез директно интегриране получваме

(3)
\begin{align} \underbrace{\int\limits_a^b f(x) dx}_{I(f)} = \underbrace{\int\limits_a^b L_{n-1}(f;x) dx}_{Q(f)} + \underbrace{\int\limits_a^b f[x_1, x_2, \dots, x_n, x] . \omega(x) dx}_{R(Q;f)} \end{align}

Тук, отново напомням, $\omega(x) = (x - x_1)(x - x_2) \dots (x - x_n)$

Така получваме и представянето на грешката (остатъка)

(4)
\begin{align} R(Q;f) = \int\limits_a^b f[x_1, x_2, \dots, x_n, x] \omega(x) dx \end{align}

Два специални случая

1.
Ако $f \in C^n[a,b]$ (т.е $f$ е достатъчно гладка) и $\omega(x)$ не си сменя знака в $[a,b]$, тогава

(5)
\begin{eqnarray} R(Q;f) & = & \int\limits_a^b f[x_1, x_2, \dots, x_n, x] \omega(x) dx = \\ & = & \int\limits_a^b \frac{f^{(n)} (\xi(x))}{n!} \omega(x) dx = \\ & = & \frac{f^{(n)} (\xi)}{n!} \int\limits_a^b \omega(x) dx \end{eqnarray}
1

2.
Ако $f \in C^n[a,b]$, $\omega(x)$ си сменя знака само в една точка от $(a,b)$ и $\int\limits_a^b \omega(x) dx = 0$
Тази точка е един от възлите, но нека я обозначим с $x_{n+1}$.
Тогава $(x - x_{n+1})\omega(x)$ не си сменя знака в $(a,b)$.

Знаем, че
$\displaystyle f[x_1, x_2, \dots, x_n, x_{n+1}, x] = \frac{f[x_1, x_2, \dots, x_n, x_{n+1}] - f[x_1, x_2, \dots, x_n, x]}{x_{n+1} - x}$
от където следва, че
$f[x_1, \dots, x_n,x] = f[x_1, \dots, x_n, x_{n+1}] + f[x_1, x_2, \dots, x_n, x_{n+1}, x](-x_{n+1} + x)$

Тогава

(6)
\begin{eqnarray} R(Q;f) & = & \int\limits_a^b \Big(\underbrace{f[x_1, x_2, \dots, x_n, x_{n+1}]}_{\mathrm{const}} + f[x_1, x_2, \dots, x_n, x_{n+1}, x](x - x_{n+1})\Big) \omega(x) dx = \\ & = & f[x_1, x_2, \dots, x_n, x_{n+1}] \underbrace{\int\limits_a^b \omega(x) dx}_{0} + \int\limits_a^b \underbrace{f[x_1, x_2, \dots, x_n, x_{n+1}, x]}_{\frac{f^{(n+1)}(\xi(x))}{(n+1)!}}(x - x_{n+1})\omega(x) dx = \\ & = & \frac{f^{(n+1)}(\xi)}{(n+1)!} \int\limits_a^b \omega(x)(x - x_{n+1}) dx \end{eqnarray}

Елементарни квадратурни формули

Квадратурна формула на правоъгълниците

Нека $n = 1$ и $x_1 = \frac{a+b}{2}$.
$\displaystyle L_0(f;x) = f\big(\frac{a+b}{2}\big) \Longrightarrow Q(f) = \int\limits_a^b L_0(f;x) dx = (b-a)f\big(\frac{a+b}{2}\big)$

Формулата
$\displaystyle Q^{\mathrm{pr}}(f) = (b - a) f \Big(\frac{a+b}{2}\Big)$
се нарича квадратурна формула на правоъгълниците

$\omega(x) = x - \frac{a+b}{2}$ - сменя си знака в $x = \frac{a+b}{2}$
Освен това,
$\displaystyle \int\limits_a^b \Big(x - \frac{a+b}{2}\Big) dx = \frac{\Big(x - \frac{a+b}{2} \Big)^2}{2} \Bigg|_a^b = 0$
Т.е попадаме във втория специален случай.
Следователно, ако $f \in C^2[a,b]$, то тогава
$\displaystyle R(Q^{\mathrm{pr}};f) = \frac{f''(\xi)}{2!} \int\limits_a^b \Big(x - \frac{a+b}{2}\Big)^2 dx$

И след като си сметнем интеграла (надявам се всеки да може)

$\displaystyle R(Q^{\mathrm{pr}}; f) = \frac{f''(\xi)}{24} (b - a)^3$

Квадратурна формула на трапеците

Нека $n = 2,\ x_1 = a$ и $x_2 = b$.
$\displaystyle L_1(f;x) = \frac{b - x}{b - a}f(a) + \frac{x - a}{b - a}f(b) \Longrightarrow Q(f) = \int\limits_a^b L_1(f;x) dx = \frac{b-a}{2}\big(f(a) + f(b)\big)$

Формулата
$\displaystyle Q^{\mathrm{tr}}(f) = \frac{b-a}{2} \big(f(a) + f(b)\big)$
се нарича квадратурна формула на трапеците

$\omega(x) = (x - a)(x - b)$ - не си сменя знака в $(a,b) \Longrightarrow$ попадаме в първия специален случай.
Следователно
$\displaystyle R(Q^{\mathrm{tr}}; f) = - \frac{f''(\xi)}{2!} \int\limits_a^b (x - a)(x - b) dx$
Т.е

$\displaystyle R(Q^{\mathrm{tr}}; f) = - \frac{f''(\xi)}{12}(b-a)^3$

NB: Ако $f''(x) > 0$ в $[a,b]$, тогава $R(Q^{\mathrm{pr}};f) > 0$ и $R(Q^{\mathrm{tr}};f) < 0$.
Изпълнено е също, че $Q^{\mathrm{pr}}(f) < I(f) < Q^{\mathrm{tr}}(f)$

Ако $f''(x) < 0$ в $[a,b]$ е изпълнено обратното.2

Квадратурна формула на Симпсон

Нека $n = 3,\ x_1 = a,\ x_2 = \frac{a+b}{2}$ и $x_3 = b$

Използваме квадратурната формула на правоъгълниците с възел $x_2$ и получаваме
$\displaystyle \int\limits_a^b f(x) dx = (b - a) f\big(\frac{a+b}{2}\big) + \frac{f''(\xi)}{24} (b - a)^3$

Сега използваме квадратурната формула на трапеците с възли $x_1,\ x_3$ и получваме
$\displaystyle \int\limits_a^b f(x) dx = \frac{b-a}{2} f(a) + \frac{b-a}{2} f(b) + \frac{f''(\eta)}{12} (b-a)^3$

Умножаваме горното уравнение с някаква константа $\alpha$, а долното - с $1 - \alpha$ и ги събираме.
$\alpha$ ще я изберем така, че квадратурната формула , която получим да е точна за $f(x) = x^2$.

Получаваме

(7)
\begin{align} \int\limits_a^b f(x) dx = \alpha(b-a) f\Big(\frac{a+b}{2}\Big) + (1 - \alpha) \frac{b - a}{2}\big( f(a) + f(b) \big) + \underbrace{\frac{(b-a)^3}{24} \big( \alpha f''(\xi) - 2(1 - \alpha)f''(\eta)\big)}_{\star} \end{align}

Искаме така получения остатък $\star$ да е 0.
Освен това, когато $f(x) = x^2$, то $f''(xi) = f''(\eta) = 2$

Следователно $\alpha$ трябва да увлетворява следното уравнение
$\displaystyle \alpha - 2(1 - \alpha) = 0 \iff \alpha = \frac{2}{3}$

Така получихме следвата квадратурна формула:
$\displaystyle Q(f) = \frac{2}{3} (b-a) f\big(\frac{a+b}{2}\big) + \frac{1}{3} \frac{b - a}{2} \big( f(a) + f(b) \big)$

След малко преобразуване…

Формулата
$\displaystyle Q^{\mathrm{s}}(f) = \frac{b-a}{6} \Big( f(a) + 4 f\big(\frac{a+b}{2}\big) + f(b) \Big)$
се нарича квадратурна формула на Симпсон

NB:
По построение $I(f) = Q^{\mathrm{s}}(f)$ за функциите $1,\ x,\ x^2$
Т.е квадратурната формула на Симпсон е с три възела и е точна за всяка функция от $\Pi_2$.
Следователно, $Q^{\mathrm{s}}(f)$ е интерполационна.

$\omega(x) = (x - a)\big(x - \frac{a+b}{2}\big)(x - b)$ - сменя си знака само в $\frac{a+b}{2}$
Освен това,
$\displaystyle \int\limits_a^b \omega(x) dx = 0 \Longrightarrow R(Q^{\mathrm{s}};f) = \frac{f^{4}(\xi)}{4!} \int\limits_a^b (x - a)\Big(x - \frac{a+b}{2}\Big)^2 (x - b) dx$
Т.е

$\displaystyle R(Q^{\mathrm{s}};f) = - \frac{f^4(\xi)}{2880} (b-a)^5$

Формула на Нютон-Коутс

Бележка:
Квадратурните формули на правоъгълниците, трапеците и Сипсон се наричат елементарни квадратурни формули.

Съставни квадратурни формули

… на правоъгълниците

Нека $m \in \mathbb{N}$. Делим интервала $[a,b]$ с точките $x_k= a + \frac{b - a}{m} k,\ k = 0, 1, \dots, m$
Т.е $x_0 = a$ и $x_m = b$.
Тогава

(8)
\begin{eqnarray} \int\limits_a^b f(x) dx & = & \sum_{k=0}^{m-1} \int\limits_{x_k}^{x_{k+1}} f(x) dx = \\ & = & \sum_{k=0}^{m-1} \bigg( \underbrace{(x_{k+1} - x_k)}_{\frac{b-a}{m}} f\Big( \frac{x_{k+1} - x_k}{2}\Big) + \frac{f''(\xi_k)}{24} (x_{k+1} - x_k)^3 \bigg) = \\ & = & \frac{b-a}{m} \sum_{k=0}^{m-1} f\big(\frac{x_{k+1} - x_k}{2}\big) + \Big(\frac{b -a}{m}\Big)^3 \frac{1}{24} \sum_{k=0}^{m-1} f''(\xi_k) \end{eqnarray}

Формулата
$\displaystyle Q^{\mathrm{pr}}_m (f) = \frac{b-a}{m} \sum_{k=0}^{m-1} f \Big(\frac{x_{k+1} + x_k}{2}\Big)$
се нарича m-та съставна квадратурна формула на правоъгълниците

Остатъкът го пооформяме малко

(9)
\begin{eqnarray} R(Q^{\mathrm{pr}}_m ; f) & = & \Big(\frac{b-a}{m}\Big)^3 \frac{1}{24} \sum_{k=0}^{m-1} f''(\xi_k) = \\ & = & \frac{(b-a)^3}{24m^2}. \underbrace{\frac{f''(\xi_0) + f''(\xi_1) + \dots + f''(\xi_{m-1}}{m}}_{f''(\zeta)} = \\ & = & \frac{f''(\zeta)}{24m^2}(b-a)^3 \end{eqnarray}

Окончателно

$\displaystyle R(Q^{\mathrm{pr}}_m ; f) = \frac{f''(\zeta)}{24m^2}(b-a)^3$

Ако $f''(\zeta) \leq M$ в $[a,b]$, то тогава $\displaystyle |R(Q^{pr}_m; f)| \leq \frac{M (b-a)^3}{24m^2}$, което клони към 0 както $\displaystyle \frac{1}{m^2}$

…на трапеците

Нека $m \in \mathbb{N}$. Разделяме интервала $[a,b]$ с точките $x_k = a + k. \frac{b-a}{m},\ k = 0,1 \dots, m$

Тогава

(10)
\begin{eqnarray} \int\limits_a^b f(x) dx & = & \sum_{k=0}^{m-1} \int\limits_{x_k}^{x_{k+1}} f(x) dx = \\ & = & \sum_{k=0}^{m-1} \bigg( \frac{x_{k+1} - x_k}{2} \Big( f(x_k) + f(x_{k+1}) \Big) - (x_{k+1} - x_k)^3 \frac{f''(\xi_k)}{12} \bigg) = \\ & = & \frac{b-a}{2m} \bigg( f(x_0) + f(x_m) + 2 \Big( f(x_1) + f(x_2) + \dots + f(x_{m-1}) \Big) \bigg) - \frac{(b-a)^3}{12m^2} \underbrace{\frac{f''(\xi_0) + f''(\xi_1) + \dots + f''(\xi_{m-1}}{m}}_{f''(\zeta)} \end{eqnarray}

Формулата
$\displaystyle Q^{\mathrm{tr}}_m (f) = \frac{b-a}{2m} \Big( f(x_0) + f(x_m) + 2 \sum_{k=1}^{m-1} f(x_k) \Big)$
се нарича m-та съставна квадратурна формула на трапеците

$\displaystyle R(Q^{\mathrm{tr}}_m; f) = - \frac{f''(\zeta)}{12m^2} (b-a)^3$

… на Симпсон

Нека $m \in \mathbb{N}$. Делим интервала $[a,b]$ с точките $x_k = a + \frac{b-a}{2m} . k,\ k = 0, 1, \dots, 2m$ (получават се 2m интервала).

Тогава

(11)
\begin{eqnarray} \int\limits_a^b f(x) dx & = & \sum_{k=0}^{m-1} \int\limits_{x_{2k}}^{x_{2k + 2}} f(x) dx = \\ & = & \sum_{k=0}^{m-1} \bigg( \frac{b-a}{6m} \Big( f(x_{2k}) + 4 f(x_{2k+1}) + f(x_{2k+2}) \Big) - \frac{f^{(4)}(\xi_k)}{2880} \Big(\frac{b-a}{m}\Big)^5 \bigg) = \\ & = & \frac{b-a}{6m} \Big( \underbrace{f(x_0)}_{a} + \underbrace{f(x_{2m})}_{b} + 2 \sum_{k=1}^{m-1} f(x_{2k}) + 4 \sum_{k=0}^{m-1} f(x_{2k+1}) \Big) - \frac{f^{(4)}(\zeta)}{2880 m^4} (b-a)^5 \end{eqnarray}

И така

Формулата
$\displaystyle Q^{\mathrm{s}}_m (f) = \frac{b-a}{6m} \Big( f(x_0) + f(x_{2m}) + 2 \sum_{k=1}^{m-1} f(x_{2k}) + 4 \sum_{k=0}^{m-1} f(x_{2k+1}) \Big)$
се нарича m-та съставна квадратурна формула на Симпсон
$\displaystyle R(Q^{\mathrm{s}}_m; f) = - \frac{f^{(4)}(\zeta)}{2880 m^4} (b-a)^5$

NB:
Ако $f \in C^4[a,b]$ и $| f^{(4)}(x) | \leq M, \ \forall x \in [a,b]$, то тогава
$\displaystyle | I(f) - Q^{\mathrm{s}}_m (f) | \leq \frac{M(b-a)^5}{2880m^4}$, което клони към 0 както $\frac{1}{m^4}$.

Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-NonCommercial-ShareAlike 3.0 License