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

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

Содержание

1 Введение
- 1.1 Постановка математической задачи
- 1.2 Построение квадратурных формул
2 Изложение метода
- 2.1 Применение квадратурных формул
- 2.2 Примеры квадратурных формул
3 Анализ метода
- 3.1 Погрешность квадратурной формулы
- 3.2 Численная устойчивость квадратурных формул
4 Числовой эксперимент
- 4.1 Исходный код функции
- 4.2 Пример
5 Рекомендации программисту
- 5.1 Автоматический выбор шага интегрирования с помощью апостериорной оценки погрешности методом Рунге
6 Заключение
7 Список литературы
8 См. также

Введение

В данной статье описывается метод вычисления собственного интеграла гладкой функции при помощи квадратурных формул. Формулы Ньютона-Котеса обладают следующими особенностями:

Необходимым условием сходимости данного метода является существование и ограниченность производной функции (порядок производной зависит от выбранной формулы)
Формулы Ньютона-Котеса обладают высоким порядком точности
Формулы порядка $n<9$ являются численно устойчивыми

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

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

( 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}},$

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

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

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

для 2 точек(метод трапеций), $m=1:$

$J_1[f]=\frac{1}{2}h(f_0+f_1),$

$D[f]=\frac{1}{12} h^3 f^{(2)}(\xi),$

для 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),$

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

Погрешность квадратурной формулы

Пусть функция $f$ имеет $(n+1)$ -ю непрерывную производную на отрезке $0 \le s \le 1$ , то есть $f(s) \in C^{(n+1)}[0,1]$ , $s_k$ - точки, в которых задано значение функции. Пусть используется квадратурная формула $n$ -го порядка. Введем функцию

$K_n(\xi)= \begin{cases} \frac {\xi ^n} {n!} ,& \xi \ge 0;\\ 0 ,& \xi <0. \end{cases}$

Тогда формула для погрешности имеет следующий вид

$D[f]=\int_0^1{F_{n+1}(t)f^{(n+1)}(t)dt},$

где

$F_{n+1}(t)=\sum_{k=0}^m {c_k K_n(s_k-t)} - \frac {(1-t)^{n+1}} {(n+1)!}$

Отсюда следует оценка для погрешности

( 8 )

$|D[f]| \le M_{n+1} z_{n+1}$

при $|f^{(n+1)}(t)| \le M_{n+1}$ , где $M_{n+1}>0$ - постоянная, и при

$z_{n+1}=\int_0^1{|F_{n+1}(t)|dt}.$

Если $F_{n+1}(t)$ не меняет знака на отрезке $0 \le s \le 1$ , то в силу теоремы о среднем имеем

( 9 )

$D[f]=f^{(n+1)}(\xi) \int_0^1 {F_{n+1}(t)dt}, \xi \in [0,1].$

Численная устойчивость квадратурных формул

Отметим, что формулы Ньютона-Котеса с $n \ge 10$ редко используются из-за их численной неустойчивости, приводящей к резкому возрастанию вычислительной погрешности. Причиной такой неустойчивости является то, что коэффициенты формул Ньютона-Котеса при больших $n$ имеют различные знаки, а именно при $n \ge 10$ существуют как положительные, так и отрицательные коэффициенты.

Рассмотрим квадратурную сумму

( 10 )

$I_n= \sum_{k=0}^n {c_k f(x_k)}.$

Предположим, что значения функции $f(x)$ , заданной на отрезке $[a,b]$ , вычисляются с некоторой погрешностью, то есть вместо точного значения получаем приближенное значение $\bar f(x_k)=f(x_k)+ \delta _k$ . Тогда вместо $I_n$ получим сумму

$\bar I_n= \sum_{k=0}^n {c_k(f(x_k)+\delta _k)}=I_n + \delta I_n,$

где

( 11 )

$\delta I_n=\sum_{k=0}^n {c_k \delta _k}.$

Поскольку квадратурная формула точна для $f(x) \equiv 1$ , имеем

( 12 )

$\sum_{k=0}^n{c_k}=\int_a^b{dx}=b-a$

и не зависит от $n$ .

Предположим теперь, что все коэффициенты $c_k$ неотрицательны. Тогда из ( 11 ) и ( 12 ) получим оценку

( 13 )

$|\delta I_n| \le \sum_{k=0}^n {|c_k||\delta _k|}=\sum_{k=0}^n{c_k|\delta _k|} \le (\mathrm{}\max_{0 \le k \le n} |\delta _k|)(b-a),$

которая означает, что при больших $n$ погрешность в вычислении квадратурной суммы ( 10 ) имеет тот же порядок, что и погрешность в вычислении функции. В этом случае говорят, что сумма ( 10 ) вычисляется устойчиво.

Если коэффициенты $c_k$ имеют различные знаки, то может оказаться, что сумма $\sum_{k=0}^n{|c_k|}$ не является равномерно ограниченной по $n$ и, следовательно погрешность в вычислении $I_n$ неограниченно возрастает с ростом $n$ . В этом случае вычисления по формуле ( 10 ) будут неустойчивы и пользоваться такой формулой при больших $n$ нельзя.

Числовой эксперимент

Приведем примеры вычисления интегралов с использованием формул Ньютона-Котеса. При реализации использовался язык C++, ниже приведен код функции, возвращающей приближенное значение интеграла.

Исходный код функции

На вход функция принимает 4 параметра

double a - левый конец исследуемого отрезка

double b - правый конец исследуемого отрезка

int Degree - степень используемого полинома

int Ndivisions - количество отрезков, на которые разбивается исходный. Совпадает со значением $k$ в формуле ( 7 )

f - интегрируемая функция

double f(double x)
{
	return 1/x;
}

double NewtonCotes(double a, double b, int Degree, int Ndivisions)
{
	int koef[10][10]={	1,0,0,0,0,0,0,0,0,0,
						1,1,0,0,0,0,0,0,0,0,
						1,4,1,0,0,0,0,0,0,0,
						1,3,3,1,0,0,0,0,0,0,
						7,32,12,32,7,0,0,0,0,0,
						19,75,50,50,75,19,0,0,0,0,
						41,216,27,272,27,216,41,0,0,0,
						751,3577,1323,2989,2989,1323,3577,751,0,0,
						989,5888,-928,10496,-4540,10496,-928,5888,989,0,
						2857,15741,1080,19344,5778,5778,19344,1080,15741,2857
						};
	double mltp[10]={1,1.0/2,1.0/3,3.0/8,2.0/45,5.0/288,1.0/140,7.0/17280,4.0/14175,9.0/89600};

	if ((Degree<0) || (Degree>9))throw "Wrong degree";
	if (a>=b) throw "Wrong segment";
	if (Ndivisions<1) Ndivisions = 1;

	double Sum,PartSum;
	double h=(b-a)/(Degree*Ndivisions);

	Sum=0;
	for (int j=0; j<Ndivisions; j++)
	{
		PartSum=0;
		for (int i=0; i<=Degree; i++)
			PartSum+= koef[Degree][i]*f(a + (i+j*Degree) * h );
		Sum+= mltp[Degree]*PartSum*h;
	}

	return Sum;
}

Пример

Вычислим по формулам Котеса степеней от 1 до 9 значение интеграла

$I=\int_1^2{ \frac{dx}{x} } = ln(2)$

Вычисления будем производить не разбивая отрезок на частичные, то есть используя $k+1$ точку, где $k$ -степень полинома. Результаты приведены ниже в таблице. Округления проводились с точностью 6 знаков после запятой.

Степень полинома	Количество узлов	Приближенное значение интеграла	Точное значение интеграла	Погрешность результата	Погрешность формулы
1	2	0.75	0.693147	0.056853	0.166667
2	3	0.694444	0.693147	0.001297	0.008333
3	4	0.69375	0.693147	0.000603	0.003704
4	5	0.693175	0.693147	0.000028	0.000372
5	6	0.693163	0.693147	0.000016	0.000210
6	7	0.693148	0.693147	0.000001	0.000026
7	8	0.693148	0.693147	0.000001	0.000016
8	9	0.693147	0.693147	0	0.000002
9	10	0.693147	0.693147	0	0

Заключение

Формулы Симпсона и Ньютона-Котеса являются хорошим аппаратом для вычисления определенного интеграла достаточное число раз непрерывно дифференцируемой функции. На примере формул $n$ -го порядка мы видим, что при ограниченной производной $(n+1)$ -го порядка точность формул есть $O(h^p)$ , где $h$ - шаг сетки, а $p>n$ . То есть при выполнении условий применимости данного метода, формулы дают высокий порядок сходимости. Однако при больших $n$ , в частности уже при $n \ge 10$ , вычисление приближенного значения интеграла становится численно неустойчиво, что делает их неприменимыми.

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

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

См. также

Источник — «http://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»

Категории: Численное интегрирование | Учебные задачи