Занятие № 14. Численное интегрирование. Квадратурные формулы Ньютона-Котеса
Задача численного интегрирования. Квадратурные формулы прямоугольников
Пусть требуется найти значение I интеграла Римана для некоторой заданной на отрезке [а, b] функции f(x). Хорошо известно, что для функций, допускающих на промежутке [а, b] конечное число точек разрыва первого рода, такое значение существует, единственно и может быть формально получено по определению:
(1)
где — произвольная упорядоченная система точек отрезка [а, b] такая, что
а ξi — произвольная точка элементарного промежутка [хi-1,_ хi].
В математическом анализе обосновывается аналитический способ нахождения значения I с помощью знаменитой формулы Ньютона-Лейбница
I = F(b)-F(a), (2)
где F(x) — некоторая первообразная для данной функции f(x).
К сожалению, применение этого весьма привлекательного подхода к вычислению I наталкивается на несколько серьезных препятствий. Самое главное из них — это несуществование первообразной среди элементарных функций для большинства элементарных функций f(x); например, таким способом не удается вычислить
Если первообразная F(x) для заданной функции f(x) все же найдена, то вычисление двух ее значений F(a) и F(b) может оказаться более трудоемким, чем вычисление существенно большего количества значений f(x).
Поскольку в общем случае значения функций находятся лишь приближенно, использование точной формулы (2) приводит к приближенному результату, который может быть более эффективно получен с помощью какой-либо специальной приближенной формулы на основе значений подынтегральной функции f(x). Такие специальные приближенные формулы для вычисления определенных интегралов называют квадратурными формулами (механическими квадратурами) или формулами численного интегрирования. Первый из этих терминов можно связать с геометрическим смыслом определенного интеграла: вычисление при f(х)>0 равносильно построению квадрата, равновеликого криволинейной трапеции с основанием [а,b] и «крышей» f(x).
Простые квадратурные формулы можно вывести непосредственно из определения интеграла, т.е. из представления (1). Зафиксировав там некоторое п≥ 1, будем иметь
(3)
Это приближенное равенство назовем общей формулой прямоугольников (площадь криволинейной трапеции приближенно заменяется площадью ступенчатой фигуры, составленной из прямоугольников, основаниями которых служат отрезки[хi-1, xi], а высотами — ординаты f(ξi), рис. 1).
Рис. 1. Геометрическая интерпретация общей формулы прямоугольников (.3)
Чтобы из общей формулы (3) получить конструктивное правило приближенного вычисления интеграла, воспользуемся свободой расположения точек х,, разбивающих промежуток интегрирования [а, b] на элементарные отрезки [xi-1, xi], и свободой выбора точек ξi на этих отрезках.
Условимся в дальнейшем пользоваться равномерным разбиением отрезка [а, b] на п частей точками хi с шагом полагая
(4)
При таком разбиении формула (3) приобретает вид
(5)
Теперь дело за фиксированием точек ξi на элементарных отрезах [xi-1, xi]. Рассмотрим три случая.
1) Положим ξi = xi-1. Тогда из (5) получаем
(6)
2) Пусть в (5) ξi = xi. Тогда имеем
(7)
Формулы (6) и (7) называются квадратурными формулами левых и правых прямоугольников соответственно. Совершенно очевидно (рис. 2), что дают двусторонние приближения к значению I интеграла от монотонной функции.
Рис. .2. Геометрическое оценивание интеграла от монотонной функции с помощью:
Можно рассчитывать на большую точность получения значения интеграла, если взять точку ξi посередине между точками xi-1 и xi. Отсюда приходим к следующему случаю.
3) Фиксируем В результате имеем квадратурную формулу средних прямоугольников или, иначе, формулу средней точки
(8)
Учитывая априорно большую точность формулы (8) по сравнению с формулами (6) и (7) при том же объёме и характере вычислений, эту симметричную формулу будем впредь называть просто формулой прямоугольников.
Полученные из определения интеграла квадратурные правила (6)—(8) не содержат в себе сведений, позволяющих сказать, каким нужно взять число п элементарных промежутков [xi-1, xi] или, иначе, каким должен быть шаг h разбиения (4) отрезка интегрирования [а, b], чтобы гарантированно найти значение I интеграла с наперед заданной точностью ε. Ответ на этот вопрос дает формула остаточного члена (глобальной погрешности) квадратурной формулы прямоугольников:
(9)
Как видно из формулы (9), при увеличении числа п элементарных отрезков, на которые разбивается промежуток интегрирования [а, b), ошибка численного интегрирования по формуле средней точки (8) убывает пропорционально квадрату шага h. Погрешность численного интегрирования непрерывно дифференцируемой функции по формулам левых и правых прямоугольников (6), (.7) убывает лишь по линейному закону.
Семейство квадратурных формул Ньютона-Котеса
Подстановка в интеграл вместо функции f(x) её интерполяционного многочлена Лагранжа той или иной степени п приводит к семейству квадратурных формул, называемых формулами Ньютона-Котеса.
Функция f(x)может быть единственным образом представлена в виде
где — интерполяционный многочлен Лагранжа,
— остаточный член,
Если система узлов интерполирования
совпадает с точками разбиения (4) отрезка [а, b] с шагом h, то замена переменной x=x0+qh трансформирует многочлен Лагранжа к виду
(10)
Для того, чтобы использовать такое выражение Ln(x) вместо f(х) в нужно изменить границы интегрирования (значению х=а соответствует значение q= 0,ах=b — значение q = п) и учесть, что dx=hdq. Таким образом, получаем
Это равенство, переписанное в виде
(11)
и есть квадратурная формула Ньютона-Котеса, где
(12)
— коэффициенты Котеса.
На самом деле, формулы (11)-(12) определяют семейство квадратурных формул. Параметром этого семейства является число п — степень интерполяционного многочлена, которым заменяется подынтегральная функция.
Для практического применения целесообразно рассмотреть некоторые частные случаи.
Пусть n = 1. Это значит, что диапазон интегрирования [a, b] содержит только два узла: x0 = a и x1 = b.
Вычислим коэффициенты Котеса для этого случая. Подставляя в (12) n = 1 и i = 0, 1 получим:
(13)
Подстановка коэффициентов (13) в (12) дает квадратурную формулу Ньютона-Котеса первого порядка:
(14)
В данном случае величина шага h равна длине всего диапазона интегрирования [a, b].
При таком приближении подынтегральная функция f (x) заменяется линейной (полиномом Лагранжа 1-й степени). Графически значение определенного интеграла I(f, a, b) изображается площадью криволинейной трапеции, ограниченной сверху (или снизу) графиком функции f (x) (см. пример на рис. 3).
Рис.3. К выводу формулы трапеций.
Сплошная линия - график интегрируемой функции f (x),
штрихпунктирная - график линейного приближения интегрируемой функции на интервале [x0, х1].
Крестиками заполнена область R, площадь которой равна погрешности квадратурной формулы (14).
Приближенное значение интеграла, выраженное квадратурной суммой (14), равно площади обычной трапеции с той же высотой h. Площадь криволинейного сектора дает погрешность R вычисления определенного интеграла с помощью квадратурной формулы (14).
Когда данный способ вычисления определенного интеграла используют на практике, то обычно диапазон интегрирования [a, b] разбивают на N одинаковых поддиапазонов шириной
h = (b - a) /N (15)
и к каждому из них применяют формулу (14). Тогда искомый интеграл представится суммой, каждый член которой будет содержать общий множитель (15). Каждое значение интегрируемой функции yk (k = 0, 1, …, N) в эту сумму будет входить дважды, кроме двух крайних y0 и yN. В результате интеграл I(f, a, b) выразится следующей суммой:
(16)
Последняя формула называется формулой трапеций.
Теперь пусть n=2. Это значит, что диапазон интегрирования [a, b] содержит три узла: х0=a, x1=(a+b)/2, x2=b. В данном случае шаг интегрирования равен половине диапазона интегрирования h=(b-a)/ 2
Коэффициенты Котеса в данном случае получаются вычислением интегралов (11) с подстановкой n = 2 для i = 0, 1, 2:
(7)
Искомый интеграл представится в виде:
(18)
где (b - a) = 2h.
Формула (18) была получена заменой подынтегральной функции f(x) на квадратичную, проходящую через заданные точки (xi, yi) , i = 0, 1, 2. Полученное значение искомого интеграла графически изображается площадью под параболой (см. рис. 4).
Рис. 4. К выводу формулы Симпсона.
Жирная сплошная линия - график интегрируемой функции f (x),
штриховая - график квадратичного приближения L2(x) интегрируемой функции на интервале [x0, x2]. Заштрихованная область R выражает погрешность квадратурной формулы (18).
Видно, что на интервале [x0, x1] квадратичное приближение дает заниженный результат, на интервале [x1, x2] - завышенный.
На практике весь диапазон интегрирования [a, b] обычно предварительно разбивается на четное число N = 2M равных интервалов, которые попарно объединяются в M поддиапазонов. Каждый поддиапазон содержит три узла: два крайних и один внутренний. К каждому из M поддиапазонов применяется формула (18), т.е. каждой тройке узлов соответствует своя интерполирующая парабола. Параболы стыкуются в совпадающих крайних узлах отдельных поддиапазонов. Номера всех таких узлов четные от i = 2 до i = 2 M - 2, поэтому при суммировании выражений (18) по поддиапазонам следует учесть, что слагаемые с этими номерами встретятся дважды.
В результате суммирования с учетом значений коэффициентов Котеса (17) получается рабочая формула
(19)
где для краткости введены обозначения следующих сумм:
(20)
Заметим, что сумма So содержит M слагаемых, а сумма Se состоит из M-1 слагаемых. Шаг интегрирования равен h=(b - a)/(2M)=(b - a)/N.
Формула (19) называется квадратурной формулой Симпсона.
При n = 3 аналогичным путем получается квадратурная формула с тремя узлами на интервале интегрирования:
(21)
где h = (b - a) / 3.
При решении прикладных задач диапазон интегрирования [a, b] предварительно разбивается на равные поддиапазоны, число которых кратно трем (N = 3M), и к каждому из поддиапазонов применяется формула (20). Тогда получится рабочая формула вида:
(22)
Шаг интегрирования равен h=(b-a)/(3M) = (b-a)/N, а отдельные суммы выражены следующими формулами:
(23)
Сумма (22) называется квадратурной формулой Ньютона или формулой «трех восьмых». В практике квадратурная формула Ньютона употребляется гораздо реже, чем формулы Симпсона и трапеций.
Задания: выполнить задание 3 ИДЗ№3.