Полиномиальная аппроксимация

Заменяя функцию f (x) некоторым обобщённым интерполяционным полиномом

(2.2)

где r (x) – остаточный член аппроксимации (ошибка аппроксимации), приходим к квадратурной формуле численного интегрирования

(2.3)

(2.4)

где величины x_i – называются узлами, c_i – весами, а R – погрешностью или остаточным членом формулы.

Лучше всего изучена замена f (x) алгебраическим полиномом. В этом случае обычно принимают, что r(х) º 1.

Формула трапеций

Заменим функцию f (x) на отрезке [ a, b ] полиномом Лагранжа первой степени с узлами x ₀ = a, x ₁ = b, т.е. уравнением прямой линии f (x) = k + m × x. Из условий

находим

(2.5)

Подставив (2.5) в (2.1), получим одну из простейших квадратурных формул – формулу трапеций

(2.6)

Для повышения точности численного интегрирования с использованием формулы (2.6) на отрезке [ a, b ] вводят достаточно густую сетку a = x ₀ < x ₁ < x ₂ < ¼ < x_N = b с шагом сетки h = x_i – x_{i -1} и к каждому шагу применяют формулу (2.6). Получают обобщённую формулу трапеций

(2.6,а)

Поскольку в данном случае порядок погрешности определён, то для повышения точности вычислений можно не только уменьшать величину шага h, но и использовать метод Рунге.

Метод Рунге

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

Из формулы (2.6,а) видно, что погрешность формулы трапеций имеет вид R = h ²y(h), где h - некоторая точка вблизи узла x_i. Если f '''(x) липшиц-непрерывна, то оценку погрешности нетрудно уточнить: R = h ²y(x_i) + O (h ³). Пусть в общем случае имеется некоторая приближённая формула z(x, h) для вычисления величины z (x) по значениям на равномерной сетке с шагом h, а остаточный член этой формулы имеет следующую структуру:

(2.7)

Произведём теперь расчёт по той же приближённой формуле для той же точки х, но используя равномерную сетку с другим шагом rh. Тогда получим значение z(x, rh), связанное с точным значением соотношением

(2.8)

Заметим, что O ((rh) ^{p +1}» O (h^{p +1}). Имея два расчёта на разных сетках, нетрудно определить величину погрешности. Для этого вычтем (2.7) из (2.8) и получим первую формулу Рунге:

(2.9)

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

Можно исключить найденную погрешность (2.9) из формулы (2.7) и получить результат с более высокой точностью по второй формуле Рунге:

(2.10)

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

Формула Симпсона

Вычислим интеграл по обобщённой формуле трапеций сначала на равномерной сетке с шагом h, а затем на сетке с вдвое более крупным шагом; вторая сетка получается из первой выбрасыванием узлов через один. Из вида остаточного члена (2.6,а) следует, что в формуле Рунге р = 2; для изменения шага сетки вдвое r = 2. Следовательно, проводя уточнение формулы трапеций для отрезка, содержащего узлы x ₀, x ₁, x ₂, получим формулу Симпсона

(2.11)

Обобщённая формула Симпсона на равномерной сетке и чётного числа узлов N имеет вид

(2.12)

а её остаточный член выражается соотношением

(2.13)

К самой формуле Симпсона, как это следует из вида её остаточного члена, также можно применить метод Рунге (здесь р = 4).

Формула средних

Если на отрезке [ a, b ] взять единственный узел квадратурной формулы х ₀, то функция аппроксимируется полиномом нулевой степени – константой f (x ₀). Поскольку симметрия формулы численного интегрирования приводит к повышению её точности, то выберем в качестве единственного узла середину отрезка интегрирования Приближённо заменяя площадь криволинейной трапеции площадью прямоугольника, получим формулу средних

(2.14)

Так же как и для формулы трапеций, для повышения точности вводится достаточно подробная сетка x_i и составляется обобщённая формула средних

(2.14,а)

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

Формула Эйлера

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

(2.15)

Если в таблице заданы значения только самой функции, а не её производных, то обобщённую формулу Эйлера можно применять, подставляя конечно-разностные выражения для f ₀', f_N '. Но эти выражения должны иметь второй порядок точности, чтобы соответствовать общей точности формулы, например

Получающиеся при этом формулы называются формулами Грегори.

 

Процесс Эйткена

У всех рассмотренных выше обобщённых формул на равномерных сетках ошибку можно разложить в ряд по степеням шага h. Значит, к ним ко всем применим метод Рунге. Но для его применения нужно знать, каков порядок точности исходной формулы.

Предположим, что порядок точности р существует, но неизвестен нам. В этом случае результаты численного интегрирования можно уточнить, если расчёты проведены на трёх (или более) сетках.

Чтобы упростить алгоритм расчёта, выберем три сетки с постоянным отношением шагов, т.е. с шагами h ₁ = h, h ₂ = qh, h ₃ = q ² h. Обозначим приближённое значение интеграла на k -той сетке через F_k и ограничимся главным членом погрешности; тогда можно записать

(2.16)

Это система трёх уравнений для определения трёх неизвестных F, a, p. Введём вспомогательные обозначения b = a h^p, g = q^p и преобразуем эту систему к следующему виду:

F – F ₁ = b, F – F ₂ = bg, F – F ₃ = bg². (2.17)

Перемножая крайние уравнения (2.17) и сравнивая результат с квадратом среднего уравнения, получим

b²g² = (F – F ₁)(F – F ₃) = (F – F ₂)²;

откуда легко получить уточнённое значение интеграла

(2.17)

Попарно вычитая уравнения (2.16) друг из друга, получим

F ₂ – F ₁ = b(1 - g), F ₃ – F ₂ = bg(1 - g),

или

q^p = g = (F ₃ – F ₂)/(F ₂ – F ₁).

Следовательно, эффективный порядок точности исходной формулы (2.15) равен

(2.18)

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

Рассмотрим вычисление интеграла 2/3. Мы не знаем, каков здесь эффективный порядок точности каждой из рассмотренных ранее формул численного интегрирования, так как у подынтегральной функции даже первая производная в нижнем пределе не ограничена, и все приведённые ранее априорные оценки погрешности здесь не применимы.

Составим таблицу значений функции (табл. 2.1) и вычислим интеграл по формулам трапеций и Симпсона при разных шагах (табл. 2.2).

Таблица 2.1. Значения функции

x	0,000	0,125	0,250	0,375	0,500	0,625	0,750	0,875	1,000
	0,00000	0,35355	0,50000	0,61237	0,70711	0,79057	0,86603	0,93541	1,00000

 

Таблица 2.2. Значения интеграла

h	Трапеций	Симпсон	Эйткен (трап)	Эйткен (Симпс)
1,000	0,50000	-	-	-
0,500	0,60355	0,63807	-	-
0,250	0,64328	0,65653	0,66801	-
0,125	0,65813	0,66308	0,66699	0,66668

Используя формулу (2.18), находим фактический порядок погрешности (или точности):а) для формулы трапеций р = 1,38 при шаге h = 0,25 и р = = 1,41 при шаге h = 0,125; б) для формулы Симпсона р = 1,5 при шаге h = = 0,125.

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

Задание № 3

Вычислить с точностью до 10^-4 заданный интеграл. Используя процесс Эйткена, оценить фактический порядок точности используемой квадратурной формулы (допускается использование процесса Эйткена при вычислении интеграла). Проверить, применим ли для заданной подынтегральной функции метод Рунге.

Вариант А

Вариант	Интеграл	Вариант	Интеграл

 

Вариант Б

Вариант	Интеграл	Вариант	Интеграл

 

Вариант В

Вариант	Интеграл	Вариант	Интеграл

 

Вариант Г

Вариант	Интеграл	Вариант	Интеграл

 

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

Пакет Maple любого релиза вычисляет этот интеграл с помощью специальной команды

>> int(1/sqrt(ln(1.0/x)), x=0.. 1.0);

1,772454

и т.д. Здесь возникают некоторые проблемы, о которых подробнее будет сказано ниже.