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

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

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

Содержание

Введение

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

  • Необходимым условием сходимости данного метода является существование и ограниченность производной функции (порядок производной зависит от выбранной формулы)
  • Формулы Ньютона-Котеса обладают высоким порядком точности
  • Формулы порядка 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]:

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

Рекомендации программисту

Автоматический выбор шага интегрирования с помощью апостериорной оценки погрешности методом Рунге

Величина погрешности численного интегрирования зависит как от шага сетки h, так и от гладкости подынтегральной функции f(x). В величину погрешности, кроме h входит также величина f^{(p)}(\xi), которая может сильно меняться на отрезке [a,b] и заранее неизвестна. Для уменьшения величины погрешности, можно измельчить сетку на заданном отрезке. Но при этом необходимо апостериорно оценивать погрешность. Такую оценку погрешности можно осуществить методом Рунге. Рассмотрим применение какой-то квадратурной формулы на частичном отрезке [x_{i-1},x_i]. Обозначим за I значение интеграла на всем отреpке, за I_i значение интеграла на i-ом частичном отрезке, за I_h приближенное значение интеграла на всем отрезке, полученное при помощи заданной квадратурной формулы и равномерном сетки с шагом h, а за I_{h,i} приближенное значение интеграла на i-ом частичном отрезке. Пусть данная квадратурная формула на данном частичном отрезке имеет порядок точности m, то есть

I_i-I_{h,i} \approx ch^m,

где с c - некоторая константа. Тогда

I_i-I_{ \frac{h}{2},i} \approx c( \frac{h}{2} )^m,

откуда получаем

I_i-I_{h,i} \approx 2^m(I_i-I_{ \frac{h}{2},i})
( 14 )
I_i-I_{ \frac{h}{2},i} \approx \frac {I_{ \frac{h}{2},i} - I_{h,i}} {2^m-1}.

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

I \approx I_h=\sum_{i=1}^N{I_{h,i}},

причем на всех частичных отрезках используются квадратурные формулы с одним и тем же порядком точности (или же,в частности, одна и та же формула). Проведем на каждом частичном отрезке [x_{i-1},x_i] вычисления дважды - один раз с шагом h_i, второй раз с шагом \frac{h_i}{2} и апостериорно оценим погрешность по правилу Рунге ( 14 ). Если для заданного \varepsilon >0 при i=1,2, \ldots N будут выполняться неравенства

( 15 )
|I_i-I_{ \frac{h}{2},i}| \approx \frac {|I_{ \frac{h}{2},i} - I_{h,i}|} {2^m-1} \le \frac {\varepsilon h_i} {b-a},

то получим

|I-I_{\frac{h}{2}}| \le \frac {\varepsilon}{b-a} \sum_{i=1}^N{h_i} \le \varepsilon,

то есть будет достигнута заданная точность \varepsilon.

Если же на каком-то из частичных отрезков оценка ( 15 ) не будет выполняться, то шаг на этом отрезке надо измельчить еще в два раза и снова оценить погрешность. Измельчение сетки на данном отрезке следует проводить до тех пор, пока не будет достигнута оценка вида ( 15 ). Заметим, что для некоторой функции f(x) такое измельчение может продолжаться слишком долго. Поэтому в соответствующей программе следует предусмотреть ограничение сверху на число измельчений.

Таким образом, автоматический выбор шага интегрирования приводит к тому, что интегрирование ведется с крупным шагом на участках плавного изменения функции f(x) и с мелким шагом - на участках быстрого изменения f(x). Это позволяет при заданной точности \epsilon уменьшить количество вычислений значений f(x) по сравнению с расчетом на сетке с постоянным шагом. Подчеркнем, что для нахождения сумм I_{ \frac{h}{2},i} не надо пересчитывать значения f(x) во всех узлах, достаточно вычислять f(x) только в новых узлах.

Заключение

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

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

См. также

Личные инструменты