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

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

Текущая версия

Введение

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

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

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

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

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

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

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

Для построения формулы численного интегрирования на всём отрезке достаточно построить квадратурную формулу для интеграла

на частичном отрезке и воспользоваться свойством (3).

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

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

Данную формулу при помощи замены можно привести к стандартному виду

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

Рассмотрим случай, когда узлы заданы и требуется найти веса квадратурной формулы . Будем пользоваться требованием: формула (5) должна быть точной для любого полинома степени . Для того чтобы полином степени удовлетворял данному требованию, достаточно потребовать, чтобы квадратурная формула была точной для любого одночлена степени . Учитывая, что , получаем уравнение

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

Так, полагая , имеем систему , решением которой являются веса формулы Симпсона: . Таким образом, формула Симпсона является точной для полинома второй степени. Однако, в силу симметрии, она является точной и для всех полиномов третьей степени:

так как она точна для :

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

где - интерполяционный коэффициент Лагранжа. Из равенства

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

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

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

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

Вернемся к рассмотрению интеграла (1). Как было показано выше, данный интеграл заменой сводится к интегралу на единичном отрезке , а следовательно легко обобщить формулы для приближенного вычисления интеграла на единичном отрезке на произвольный. Применим аппарат квадратурных формул. Пусть задано равномерное разбиение отрезка с шагом , где обозначим . Пусть также выбрана некоторая квадратурная формула Ньютона-Котеса (то есть выбрана степень полинома , а следовательно каждый полином строится по точке сетки). Мы также считаем, что данный набор из точки можно разбить на поднаборы по точек с совпадающими крайним, то есть

где

Тогда, суммируя значения квадратур на каждом поднаборе, получим приближенное значение искомого интеграла. Если обозначить весовые множители то приближенное значение интеграла можно записать в виде двойной суммы

Данный алгоритм естественным образом обобщается на случай, когда , где , а на каждом отрезке задана равномерная сетка. Тогда искомый интеграл равен

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

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

Приведем примеры квадратурных формул Котеса на равномерной сетке с шагом , где обозначим  :

• для 2 точек(метод трапеций),

• для 3 точек(метод Симпсона),

• для 4 точек,

• для 5 точек,

• для 6 точек,

• для 7 точек,

• для 8 точек,

• для 9 точек,

• для 10 точек,

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

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

Пусть функция имеет -ю непрерывную производную на отрезке , то есть , - точки, в которых задано значение функции. Пусть используется квадратурная формула -го порядка. Введем функцию

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

где

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

при , где - постоянная, и при

Если не меняет знака на отрезке , то в силу теоремы о среднем имеем

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

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

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

Предположим, что значения функции , заданной на отрезке , вычисляются с некоторой погрешностью, то есть вместо точного значения получаем приближенное значение . Тогда вместо получим сумму

где

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

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

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

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

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

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

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

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

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

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

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

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

int Ndivisions - количество отрезков, на которые разбивается исходный. Совпадает со значением в формуле ( 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 значение интеграла

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

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

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

Величина погрешности численного интегрирования зависит как от шага сетки , так и от гладкости подынтегральной функции . В величину погрешности, кроме входит также величина , которая может сильно меняться на отрезке и заранее неизвестна. Для уменьшения величины погрешности, можно измельчить сетку на заданном отрезке. Но при этом необходимо апостериорно оценивать погрешность. Такую оценку погрешности можно осуществить методом Рунге. Рассмотрим применение какой-то квадратурной формулы на частичном отрезке . Обозначим за значение интеграла на всем отреpке, за значение интеграла на -ом частичном отрезке, за приближенное значение интеграла на всем отрезке, полученное при помощи заданной квадратурной формулы и равномерном сетки с шагом , а за приближенное значение интеграла на -ом частичном отрезке. Пусть данная квадратурная формула на данном частичном отрезке имеет порядок точности , то есть

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

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

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

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

то получим

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

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

Таким образом, автоматический выбор шага интегрирования приводит к тому, что интегрирование ведется с крупным шагом на участках плавного изменения функции и с мелким шагом - на участках быстрого изменения . Это позволяет при заданной точности уменьшить количество вычислений значений по сравнению с расчетом на сетке с постоянным шагом. Подчеркнем, что для нахождения сумм не надо пересчитывать значения во всех узлах, достаточно вычислять только в новых узлах.

Заключение

Формулы Симпсона и Ньютона-Котеса являются хорошим аппаратом для вычисления определенного интеграла достаточное число раз непрерывно дифференцируемой функции. На примере формул -го порядка мы видим, что при ограниченной производной -го порядка точность формул есть , где - шаг сетки, а . То есть при выполнении условий применимости данного метода, формулы дают высокий порядок сходимости. Однако при больших , в частности уже при , вычисление приближенного значения интеграла становится численно неустойчиво, что делает их неприменимыми.

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

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

См. также

• Практикум ММП ВМК, 4й курс, осень 2008
• Код программы