6.Численное интегрирование При вычислении определенных интегралов с помощью формулы Ньютона- Лейбница (6.1) необходимо для подынтегральной функции f(x) найти первообразную F(x). Если интеграл не выражается через элементарные функции (тогда соответствующий неопределенный интеграл определяет новую функцию), то для его вычисления используют численные методы. Примеры таких интегралов приведены ниже: Численные методы интегрирования применяют тогда, когда аналитические методы интегрирования не применимы или слишком сложны. Если необходимо вычислить определенный интеграл от таблично заданной функции, то применение численного интегрирования является неизбежным. Формулы для приближенного вычисления интегралов часто называют квадратурными формулами.

6.1.Квадратурные формулы Ньютона-Котеса Пусть требуется вычислить определенный интеграл (6.2) При выводе квадратурных формул для приближенного вычисления определенного интеграла вспомним его геометрический смысл интеграл равен площади криволинейной трапеции, ограниченной графиком подынтегральной функции, осью Ох и отрезками прямых х = а и х = b. Разобьем отрезок [a, b] на n частей точками x i : (6.3) Обозначим через y i значения функции в точках x i. Заменим подынтегральную функцию интерполяционным многочленом Лагранжа (4.17): (6.4)

Тогда получим приближенную формулу для вычисления интеграла: (6.5) где A i числовые коэффициенты, которые не зависят от подынтегральной функции и их значения можно определить для заданного n. Выведем формулы для вычисления коэффициентов A i. Введем обозначения Тогда многочлен Лагранжа можно записать в виде Заменяя под знаком интеграла в (6.5) функцию y(x) многочленом L n (x), получим

Отсюда следуют формулы для вычисления коэффициентов A i : Коэффициенты H i называются коэффициентами Котеса и вычисляются по формулам: (6.6) В следующих пунктах рассмотрим простейшие квадратурные формулы, являющиеся частными случаями этих соотношений и получающиеся таким способом.

Формула прямоугольников Формулы прямоугольников получаются заменой подынтегральной функции постоянным значением. В качестве такого значения выбирают значение функции в одной из точек отрезка [a, b], или на левом конце отрезка, или на правом конце отрезка, или в середине отрезка (рис.6.1): (6.7) (6.8) (6.9) рис. 6.1 ababab y xxx yy а) б) в)

Если формулы (6.6) (6.8) применить к каждой части [x i, x i + 1 ] отрезка [a, b], то получим общие формулы прямоугольников. Фактически, определенный интеграл приближенно заменяется интегральной суммой: (6.10) (6.11) (6.12) Геометрически это означает, что площадь криволинейной трапеции приближенно заменяется площадью ступенчатой фигуры. В частности, рис.6.2 иллюстрирует формулу (6.10). рис. 6.2 ab x y

Формулу (6.10) называют формулой левых прямоугольников, а формулу (6.11) формулой правых прямоугольников, а (6.12), соответственно, формулой средних прямоугольников. Формулы прямоугольников практически не используются из-за большой погрешности порядка O(h) (у формулы средних более высокий прядок O(h 2 )). Положим в формулах (6.6) n = 1 и вычислим значения A i : Мы заменили подынтегральную функцию многочленом Лагранжа первой степени и получили формулу трапеций: (6.13) Формула трапеций

Геометрический смысл формулы трапеций (6.13) заключается в том, что кривая y = y(x) заменяется отрезком прямой, проходящей через точки (x 0, y 0 ) и (x 1, y 1 ), или, в других обозначениях, (a, y(a)) и (b, y(b)) (рис.6.3). Рис. 6.3 Заметим, что формулы трапеций и средних прямоугольников являются точными для линейной функции. Если обобщить (6.13) для равномерного разбиения отрезка на n частей, то приходим к общей формуле трапеций (рис.6.4): (6.14) Рис ab x y ab x y

Погрешность формулы трапеций (6.13) есть величина порядка O(h3). В этом можно убедиться, используя формулу погрешности интерполяционной формулы Лагранжа. А для общей формулы трапеций (6.14) погрешность есть величина порядка O(h 2 ), так как при суммировании погрешности накапливаются. Применяя интерполяционную формулу Лагранжа при n = 2, получим значения коэффициентов: (6.15) Геометрический смысл формулы Симпсона (6.15) заключается в том, что кривая y = y(x) заменяется частью параболы, проходящей через три точки (x 0, y 0 ), (x 1, y 1 ) и (x 2, y 2 ) Формула Симпсона

Формула Симпсона точна не только для полинома второй степени, но и для полинома третьей степени в силу симметрии, показанной на рис. 6.4б. Рис.6.4 Остаточный член формулы (6.15) имеет порядок O(h 5 ). Общая формула Симпсона строится для четного n = 2m; при этом формула (6.15) применяется для каждого из отрезков [x 0, x 2 ], [x 2, x 4 ], …, [x n–2, x n ]: (6.16) Остаточный член общей формулы Симпсона (6.16) имеет порядок O(h 4 ). ab x y ab x y - + а)б) y=y(x) y=L 2 (x) y=L 3 (x) y=L 2 (x)

Формулы Ньютона-Котеса высших порядков Полагая в формулах (6.6) n = 3, можно вычислить значения коэффициентов A i и получить квадратурную формулу Ньютона (6.17) Формула (6.17) имеет погрешность того же порядка, что и формула Симпсона (6.15), т.е. O(h 5 ). Приведем таблицу значений коэффициентов Котеса (табл.6.3). Таблица 6.3

Таблица 6.3 (продолжение) Например, квадратурная формула Ньютона-Котеса для n = 5 имеет вид: (6.18) Квадратурные формулы с нечетным числом узлов (n = 2, 4, 6) являются более удобными, т.е. выгоднее применять формулу Симпсона (6.15), чем формулу Ньютона (6.17), так как при одном и том же порядке погрешности, формула Ньютона требует больше узлов (и больше вычислений), чем формула Симпсона. Аналогично, формула для n = 4 лучше, чем формула для n = 5 и т.д.

6.2.Правило Рунге оценки погрешности Приведем сводку квадратурных формул с остаточными членами. Формула трапеций: Формула Симпсона: Рассмотрим на примере общей формулы Симпсона (6.16) правило Рунге для оценки погрешности. Пусть подынтегральная функция четырежды непрерывно дифференцируема. Тогда формула остаточного члена имеет вид: (6.19) где ξ некоторое число из отрезка [a, b]. Предположим, что производная y IV (x) изменяется на этом отрезке медленно, и приближенно можно записать остаточный член в виде где M постоянная. Пусть S h и S 2h соответственно значения интеграла, полученные по общей формуле Симпсона с шагом h и 2h. Тогда справедливы соотношения

Отсюда получим формулу для оценки погрешности (6.20) Уточненное значение интеграла по формуле Симпсона вычисляется с учетом поправки (6.21) Для формулы трапеций также можно применить правило Рунге. Так как формула остаточного члена общей формулы трапеций может быть представлена в виде, то справедливы соотношения

Отсюда получим формулу для оценки погрешности (6.22) Теперь для интеграла можно записать по правилу Рунге формулу трапеций с поправкой (6.23) Замечание. Здесь необходимо подчеркнуть, что описанное правило Рунге применимо только тогда, когда выполняются указанные выше условия для функции y(x) существование производной соответствующего порядка и её ограниченность, точнее говоря, возможность приближенного представления погрешности в виде для формулы Симпсона ( для формулы трапеций), где M постоянная. Погрешность представления остаточного члена в указанном виде здесь считается достаточно малой. Эти условия для конкретной функции могут не выполняться, тогда правило Рунге не будет гарантировать приемлемый результат.

6.3.Быстропеременные функции. Метод Филона Рассмотренные выше методы интегрирования хороши, если подынтегральная функция и её производные непрерывны и ограничены на отрезке интегрирования. В радиотехнических задачах часто встречается интеграл (6.23) где высокочастотное колебание (т.е. ω велико), а f(x) амплитуда. Подынтегральные функции в (6.23) являются быстропеременными (быстро осциллирующими), а их m-ая производная есть величина порядка O(ωm). Функция имеет на отрезке [a, b] примерно ω(b – a)/π корней. Так как число корней велико, то для приближенного вычисления интеграла с помощью квадратурных формул придется использовать интерполяционные многочлены высокого порядка или очень большое число узлов интерполирования, что приводит к громоздким вычислениям и, как следствие, большим ошибкам округления.

В качестве узлов интерполирования примем x i = (b + a)/2 + (b – a)d i /2, i = 1, …, n, (6.24) где числа d i принадлежат отрезку [–1; 1]. Заменим функцию в интеграле (6.23) на интерполяционный многочлен с узлами x i : (6.25) Здесь (6.26) Мы получили квадратурную формулу (6.27) Рассмотрим частные случаи. При n = 2, d 1 = –1, d 2 = 2 из (6.26), (6.27) получим: (6.28)

(6.29) (6.30) Проведем анализ влияния погрешности машинных округлений на значения D1(p), D2(p), вычисляемых по формулам (6.28), (6.29). Пусть t число двоичных разрядов записи числа. Тогда можно считать, что погрешность вычисления значений sinp, cosp есть величина порядка O(2 –t ), а погрешность вычисления D 1 (p), D 2 (p) имеет порядок O(2 –t /p). Отсюда следует, что при малых значениях p погрешности D1(p), D2(p) могут быть большими. С другой стороны, используя разложения в ряд Тейлора можно вычислить пределы Чтобы обойти эту ситуацию при вычислении D1(p), D2(p) в стандартных программах применяется следующее правило: (6.31)

где p 2 некоторое малое положительное число, определяемое подбором. Разделим отрезок [a, b] на N частей точками x i = a + ih, i = 0, 1, …, N; h = (b – a)/N. Тогда можем записать обобщенную формулу (6.30) в виде: (6.32) Сформулируем алгоритм вычисления интеграла (6.23) по формулам (6.31) (6.32). Алгоритм метода Филона. 1. Зададим N. Вычислим узловые точки: h = (b – a)/N; x i = a + ih, i = 0, 1, …, N; (6.33) 2. Вычислим параметры: ; p 2 = 0,00001; (6.34)

(6.35) (6.36) 3. Вычислим сумму: (6.37) Пример 6.5. Вычислить интеграл. Решение в Mathcad. Интеграл можно вычислить точно интегрированием по частям:

Вычислим значение полученного выражения в программе Mathcad (Обратите внимание, i мнимая единица, для её ввода нужно нажать указателем мыши на букву i на панели инструментов «Калькулятор»): Интересно отметить, что встроенная программа вычисления интеграла дает практически тот же результат Применим алгоритм (6.33) (6.37) в Mathcad:

Как видим, результат метода Филона дает четыре верных знака после запятой при числе узлов N = 320. Вычислим число нулей подынтегральной функции: На каждую полуволну приходится примерно два узла. Воспользуемся формулой трапеций при тех же узлах: Формула трапеций дает неверный результат. Если мы увеличим число узлов до 1600, то метод Филона дает пять верных знаков после запятой, а метод трапеций всего два. Замечание. Если амплитуда f(x) = 1, то метод Филона дает точное значение интеграла при любом числе узлов, так как в этом случае на каждом промежуточном отрезке применяется формула Ньютона-Лейбница и не используется замена подынтегральной функции приближенной.