Методы парабол (Симпсона) и более высоких степеней (Ньютона - Котеса)

Материал из MachineLearning.

(Различия между версиями)

Перейти к: навигация, поиск

Версия 19:36, 16 октября 2008

Содержание

1 Введение
- 1.1 Постановка математической задачи
- 1.2 Построение квадратурных формул
2 Изложение метода
- 2.1 Применение квадратурных формул
- 2.2 Примеры квадратурных формул
3 Анализ метода
4 Числовой пример
5 Рекомендации программисту
6 Заключение
7 Список литературы

Введение

Постановка математической задачи

Задача численного интегрирования состоит в нахождении приближенного значения интеграла

( 1 )

$J[f]=\int_a^b{f(x)dx},$

где $f(x)$ - заданная и интегрируемая на отрезке $[a,b]$ функция. На отрезке вводится сетка $\omega=\{x_i:x_0=a<x_1<\ldots<x_i<\ldots<x_N=b\}$ и в качестве приближенного значения интеграла рассматривается число

( 2 )

$J_N[f]=\sum_{i=0}^N {c_i f(x_i)},$

где $f(x_i)$ - значения функции $f(x)$ в узлах $x=x_i$ , где $c_i$ - весовые множители, зависящие только от узлов, но не зависящие от выбора $f(x)$ . Формула (2) называется квадратурной формулой. Задача численного интегрирования при помощи квадратур состоит в отыскании таких узлов $\{x_i\}$ и таких весов $\{c_i\}$ , чтобы погрешность квадратурной формулы

$D[f]=\sum_{i=0}^N{c_i f(x_i)} - \int_a^b{f(x)dx} = J_N[f] - J[f]$

была минимальной по модулю для функции из заданного класса (величина $D[f]$ зависит от гладкости $f(x)$ ). Погрешность зависит как от расположения узлов, так и от выбора весовых коэффициентов. Введем на $[a,b]$ равномерную сетку с шагом $h$ , т.е. множество точек $\omega_h=\{x_i=a+ih, i=0,1,\ldots,N,hN=b-a}$ , и представим интеграл (1) в виде суммы интегралов по частичным отрезкам:

( 3 )

$\int_a^b{f(x)dx}=\sum_{i=1}^N{\int_{x_{i-1}}^{x_i}{f(x)dx}}.$

Для построения формулы численного интегрирования на всм отрезке $[a,b]$ достаточно построить квадратурную формулу для интеграла

( 4 )

$\int_{x_{i-1}}^{x_i}{f(x)dx}$

на частичном отрезке $[x_{i-1},x_i]$ и воспользоваться свойством (3).

Построение квадратурных формул

В силу вышеизложенного выше, вычисление приближенного значения интеграла производится при помощи квадратурной формулы

$\int_a^b{f(s)ds}\approx\sum_{k=0}^m{c_kf(s_k)}.$

Данную формулу при помощи замены $x=a+(b-a)s, f(x)=f(a+(b-a)s)$ можно привести к стандартному виду

( 5 )

$\int_0^1{f(s)ds}\approx\sum_{k=0}^m{c_kf(s_k)}.$

В общем случае узлы и веса неизвестны и подлежат определению.

Рассмотрим случай, когда узлы заданы и требуется найти веса квадратурной формулы $\{c_k\}$ . Будем пользоваться требованием: формула (5) должна быть точной для любого полинома $P_r(s)$ степени $r\le m$ . Для того чтобы полином степени $r$ удовлетворял данному требованию, достаточно потребовать, чтобы квадратурная формула была точной для любого одночлена $s^\sigma$ степени $\sigma (\sigma=0,1,\ldots,r)$ . Учитывая, что $J[s^\sigma]=\frac{1}{\sigma+1}$ , получаем $m+1$ уравнение

$c_0+c_1+\ldots+c_m=1,$

$c_0 s_0+c_1 s_1+\ldots+c_m s_m=\frac{1}{2},$

$.........................................$

$c_0 s_0^\sigma + c_1 s_1^\sigma + \ldots + c_m s_m^\sigma=\frac{1}{\sigma+1},$

$.........................................$

$c_0 s_0^m + c_1 s_1^m + \ldots + c_m s_m^m=\frac{1}{m+1}.$

Эта система имеет единственное решение, так как ее определителем является определитель Вандермонда, отличный от нуля если нет совпадающих узлов, $s_0<s_q<\ldots<s_m$ .

Так, полагая $m=2,s_0=0,s_1=\frac{1}{2},s_2=1$ , имеем систему $c_0+c_1+c_2=1, \frac{c_1}{2}+c_2=\frac{1}{2}, \frac{c_1}{4}+c_2=\frac{1}{3}$ , решением которой являются веса формулы Симпсона: $c_0=c_2=\frac{1}{6},c_1=\frac{4}{6}$ . Таким образом, формула Симпсона является точной для полинома второй степени. Однако, в силу симметрии, она является точной и для всех полиномов третьей степени:

$P_3(s)=1+\alpha_1 (s-\frac{1}{2})+\alpha_2 (s-\frac{1}{2})^2 +\alpha_3 (s-\frac{1}{2})^3,$

так как она точна для $f(s)=(s-\frac{1}{2})^3$ :

$J_2[(s-\frac{1}{2})^3]=\frac{1}{6}((-\frac{1}{2})^3+4\cdot 0+(\frac{1}{2})^3)=0,$

$L[(s-\frac{1}{2})^3]=\int_0^1{(s-\frac{1}{2})^3ds}=0.$

Формулы треугольника и трапеции точны для линейной функции,т.е. для полинома первой степени, в чем легко убедиться непосредственно. В общем случае в качестве $P_m(s)$ можно выбрать интерполяционный полином Лагранжа

$P_m(s)=\sum_{k=0}^m{l_k^{(m)}(s)f(s_k)},$

где $l_k^{(m)}(s)$ - интерполяционный коэффициент Лагранжа. Из равенства

$L[P_m]=\int_0^1{P_m(s)ds}=\sum_{k=0}^m{f(s_k)\int_0^1{l_k^{(m)}(s)ds}}=\sum_{k=0}^m{c_k f(s_k)}$

видно, что формула (5) верна для полинома степени $m$ , если весовые коэффициенты $c_k$ определяются по формуле

( 6 )

$c_k=\int_0^1{l_k^{(m)}(s)ds}.$

Формулы такого типа называют квадратурными формулами Котеса.

Изложение метода

Применение квадратурных формул

Вернемся к рассмотрению интеграла (1). Как было показано выше, данный интеграл заменой сводится к интегралу на единичном отрезке $[0,1]$ , а следовательно легко обобщить формулы для приближенного вычисления интеграла на единичном отрезке на произвольный. Применим аппарат квадратурных формул. Пусть задано равномерное разбиение отрезка с шагом $h:$ $\omega_h=\{x_i=a+ih, i=0,1,\ldots,N;hN=b-a}$ , где обозначим $f_i=f(a+ih)$ . Пусть также выбрана некоторая квадратурная формула Ньютона-Котеса (то есть выбрана степень полинома $n$ , а следовательно каждый полином строится по $n+1$ точке сетки). Мы также считаем, что данный набор из $N+1$ точки можно разбить на поднаборы по $m$ точек с совпадающими крайним, то есть

( 7 )

$N=n*k,$ где $k \in \mathbb N.$

Тогда, суммируя значения квадратур на каждом поднаборе, получим приближенное значение искомого интеграла. Если обозначить весовые множители $c_i,i=0 \ldots n,$ то приближенное значение интеграла можно записать в виде двойной суммы

$\int_a^b{f(x)dx \approx \sum_{j=0}^{k-1}{\sum_{i=0}^n{c_i f_{(i+j*n)}}}$

Данный алгоритм естественным образом обобщается на случай, когда $[a,b]=\bigcup_{i=0}^K{[a_i,b_i]}$ , где $a_0=a, b_K=b, b_{i+1}=a_i;i=0\ldots K-1$ , а на каждом отрезке $[a_i,b_i]$ задана равномерная сетка. Тогда искомый интеграл равен

$\int_a^b{f(x)dx} = \sum_{i=0}^K{\int_{a_i}^{b_i}{f(x)dx}},$

а на каждом из отрезков приближенное значение интеграла вычисляется при помощи квадратурных формул.

Примеры квадратурных формул

Приведем примеры квадратурных формул Котеса на равномерной сетке с шагом $h:$ $\omega_h=\{x_i=a+ih, i=0,1,\ldots,N;hN=b-a}$ , где обозначим $f_i=f(a+ih),$ $\xi \in [a,b]$ :

для 3 точек(метод Симпсона), $m=2:$

$J_2[f]=\frac{1}{3}h(f_0+4f_1+f_2),$

$D[f]=\frac{1}{90} h^5 f^{(4)}(\xi),$

для 4 точек, $m=3:$

$J_3[f]=\frac{3}{8}h(f_0+3f_1+3f_2+f_3),$

$D[f]=\frac{3}{80} h^5 f^{(4)}(\xi),$

для 5 точек, $m=4:$

$J_4[f]=\frac{2}{45}h(7f_0+32f_1+12f_2+32f_3+7f_4),$

$D[f]=\frac{8}{945} h^7 f^{(6)}(\xi),$

для 6 точек, $m=5:$

$J_5[f]=\frac{5}{288}h(19f_0+75f_1+50f_2+50f_3+75f_4+19f_5),$

$D[f]=\frac{275}{12096} h^7 f^{(6)}(\xi),$

для 7 точек, $m=6:$

$J_6[f]=\frac{1}{140}h(41f_0+216f_1+27f_2+272f_3+27f_4+216f_5+41f_6),$

$D[f]=\frac{9}{1400} h^9 f^{(8)}(\xi),$

для 8 точек, $m=7:$

$J_7[f]=\frac{7}{17280}h(751f_0+3577f_1+1323f_2+2989f_3+2989f_4+1323f_5+3577f_6+751f_7),$

$D[f]=\frac{8183}{518400} h^9 f^{(8)}(\xi),$

для 9 точек, $m=8:$

$J_8[f]=\frac{4}{14175}h(989f_0+5888f_1-928f_2+10496f_3-4540f_4+10496f_5-928f_6+5888f_7+989f_8),$

$D[f]=\frac{2368}{467775} h^{11} f^{(10)}(\xi),$

для 10 точек, $m=9:$

$J_9[f]=\frac{9}{89600}h(2857f_0+15741f_1+1080f_2+19344f_3+5778f_4+5778f_5+19344f_6+1080f_7+15741f_8+2857f_9).$

$D[f]=\frac{173}{14620} h^{11} f^{(10)}(\xi),$

Анализ метода

Числовой пример

Заключение

Список литературы

А.А.Самарский, А.В.Гулин. Численные методы М.: Наука, 1989.
А.А.Самарский. Введение в численные методы М.: Наука, 1982.
http://mathworld.wolfram.com/Newton-CotesFormulas.html

Это незавершённая статья. Вы поможете проекту, исправив и дополнив её.

Источник — «http://www.machinelearning.ru/wiki/index.php?title=%D0%9C%D0%B5%D1%82%D0%BE%D0%B4%D1%8B_%D0%BF%D0%B0%D1%80%D0%B0%D0%B1%D0%BE%D0%BB_%28%D0%A1%D0%B8%D0%BC%D0%BF%D1%81%D0%BE%D0%BD%D0%B0%29_%D0%B8_%D0%B1%D0%BE%D0%BB%D0%B5%D0%B5_%D0%B2%D1%8B%D1%81%D0%BE%D0%BA%D0%B8%D1%85_%D1%81%D1%82%D0%B5%D0%BF%D0%B5%D0%BD%D0%B5%D0%B9_%28%D0%9D%D1%8C%D1%8E%D1%82%D0%BE%D0%BD%D0%B0_-_%D0%9A%D0%BE%D1%82%D0%B5%D1%81%D0%B0%29»

Категории: Незавершённые статьи | Численное интегрирование

@@ Строка 2: / Строка 2: @@
 === Постановка математической задачи ===
-:Задача численного интегрирования состоит в нахождении приближенного значения интеграла
+Задача численного интегрирования состоит в нахождении приближенного значения интеграла
 {{ eqno | 1 }}
@@ Строка 13: / Строка 13: @@
 где <tex>f(x_i)</tex> - значения функции <tex>f(x)</tex> в узлах <tex>x=x_i</tex> , где <tex>c_i</tex> - ''весовые множители'', зависящие только от узлов, но не зависящие от выбора <tex>f(x)</tex>. Формула {{eqref|2}} называется ''квадратурной формулой''.
-:Задача численного интегрирования при помощи квадратур состоит в отыскании таких узлов <tex>\{x_i\}</tex> и таких весов <tex>\{c_i\}</tex>, чтобы ''погрешность квадратурной формулы''
+Задача численного интегрирования при помощи квадратур состоит в отыскании таких узлов <tex>\{x_i\}</tex> и таких весов <tex>\{c_i\}</tex>, чтобы ''погрешность квадратурной формулы''
 ::<tex>D[f]=\sum_{i=0}^N{c_i f(x_i)} - \int_a^b{f(x)dx} = J_N[f] - J[f]</tex>
@@ Строка 32: / Строка 32: @@
 === Построение квадратурных формул ===
-:В силу вышеизложенного выше, вычисление приближенного значения интеграла производится при помощи квадратурной формулы
+В силу вышеизложенного выше, вычисление приближенного значения интеграла производится при помощи квадратурной формулы
 ::<tex>\int_a^b{f(s)ds}\approx\sum_{k=0}^m{c_kf(s_k)}.</tex>
@@ Строка 85: / Строка 85: @@
 Формулы такого типа называют квадратурными '''формулами Котеса'''.
+== Изложение метода ==
+=== Применение квадратурных формул ===
+Вернемся к рассмотрению интеграла {{eqref|1}}. Как было показано выше, данный интеграл заменой сводится к интегралу на единичном отрезке <tex>[0,1]</tex>, а следовательно легко обобщить формулы для приближенного вычисления интеграла на единичном отрезке на произвольный. Применим аппарат квадратурных формул. Пусть задано равномерное разбиение отрезка с шагом <tex>h:</tex> <tex>\omega_h=\{x_i=a+ih, i=0,1,\ldots,N;hN=b-a}</tex>, где обозначим <tex>f_i=f(a+ih)</tex>. Пусть также выбрана некоторая квадратурная формула Ньютона-Котеса (то есть выбрана степень полинома <tex>n</tex>, а следовательно каждый полином строится по <tex>n+1</tex> точке сетки). Мы также считаем, что данный набор из <tex>N+1</tex> точки можно разбить на поднаборы по <tex>m</tex> точек с совпадающими крайним, то есть
+{{ eqno | 7 }}
+::<tex>N=n*k,</tex> где <tex>k \in \mathbb N.</tex>
+Тогда, суммируя значения квадратур на каждом поднаборе, получим приближенное значение искомого интеграла. Если обозначить весовые множители <tex>c_i,i=0 \ldots n,</tex> то приближенное значение интеграла можно записать в виде двойной суммы
+::<tex>\int_a^b{f(x)dx \approx \sum_{j=0}^{k-1}{\sum_{i=0}^n{c_i f_{(i+j*n)}}}</tex>
+Данный алгоритм естественным образом обобщается на случай, когда <tex>[a,b]=\bigcup_{i=0}^K{[a_i,b_i]}</tex>, где <tex>a_0=a, b_K=b, b_{i+1}=a_i;i=0\ldots K-1</tex>, а на каждом отрезке <tex>[a_i,b_i]</tex> задана равномерная сетка. Тогда искомый интеграл равен
+::<tex>\int_a^b{f(x)dx} = \sum_{i=0}^K{\int_{a_i}^{b_i}{f(x)dx}},</tex>
+а на каждом из отрезков приближенное значение интеграла вычисляется при помощи квадратурных формул.
 === Примеры квадратурных формул ===
-:Приведем примеры квадратурных формул Котеса на равнемерной сетке с шагом <tex>h:</tex> <tex>\omega_h=\{x_i=a+ih, i=0,1,\ldots,N,hN=b-a}</tex>, где обозначим <tex>f_i=f(a+ih)</tex>:
+Приведем примеры квадратурных формул Котеса на равномерной сетке с шагом <tex>h:</tex> <tex>\omega_h=\{x_i=a+ih, i=0,1,\ldots,N;hN=b-a}</tex>, где обозначим <tex>f_i=f(a+ih),</tex>&nbsp;<tex>\xi \in [a,b]</tex>:
 *для 3 точек(метод Симпсона), <tex>m=2:</tex>
 ::<tex>J_2[f]=\frac{1}{3}h(f_0+4f_1+f_2),</tex>
-::<tex>c_0=c_2=\frac{1}{6},c_1=\frac{4}{6}.</tex>
+::<tex>D[f]=\frac{1}{90} h^5 f^{(4)}(\xi),</tex>
 *для 4 точек, <tex>m=3:</tex>
 ::<tex>J_3[f]=\frac{3}{8}h(f_0+3f_1+3f_2+f_3),</tex>
-::<tex>c_0=c_3=\frac{1}{8},c_1=c_2=\frac{3}{8}.</tex>
+::<tex>D[f]=\frac{3}{80} h^5 f^{(4)}(\xi),</tex>
 *для 5 точек, <tex>m=4:</tex>
 ::<tex>J_4[f]=\frac{2}{45}h(7f_0+32f_1+12f_2+32f_3+7f_4),</tex>
-::<tex>c_0=c_4=\frac{7}{90},c_1=c_3=\frac{32}{90},c_2=\frac{12}{90}.</tex>
+::<tex>D[f]=\frac{8}{945} h^7 f^{(6)}(\xi),</tex>
 *для 6 точек, <tex>m=5:</tex>
 ::<tex>J_5[f]=\frac{5}{288}h(19f_0+75f_1+50f_2+50f_3+75f_4+19f_5),</tex>
+::<tex>D[f]=\frac{275}{12096} h^7 f^{(6)}(\xi),</tex>
 *для 7 точек, <tex>m=6:</tex>
 ::<tex>J_6[f]=\frac{1}{140}h(41f_0+216f_1+27f_2+272f_3+27f_4+216f_5+41f_6),</tex>
+::<tex>D[f]=\frac{9}{1400} h^9 f^{(8)}(\xi),</tex>
 *для 8 точек, <tex>m=7:</tex>
 ::<tex>J_7[f]=\frac{7}{17280}h(751f_0+3577f_1+1323f_2+2989f_3+2989f_4+1323f_5+3577f_6+751f_7),</tex>
+::<tex>D[f]=\frac{8183}{518400} h^9 f^{(8)}(\xi),</tex>
 *для 9 точек, <tex>m=8:</tex>
 ::<tex>J_8[f]=\frac{4}{14175}h(989f_0+5888f_1-928f_2+10496f_3-4540f_4+10496f_5-928f_6+5888f_7+989f_8),</tex>
+::<tex>D[f]=\frac{2368}{467775} h^{11} f^{(10)}(\xi),</tex>
 *для 10 точек, <tex>m=9:</tex>
 ::<tex>J_9[f]=\frac{9}{89600}h(2857f_0+15741f_1+1080f_2+19344f_3+5778f_4+5778f_5+19344f_6+1080f_7+15741f_8+2857f_9).</tex>
+::<tex>D[f]=\frac{173}{14620} h^{11} f^{(10)}(\xi),</tex>
-== Изложение метода ==
 == Анализ метода ==