Численное интегрирование функций одной переменной.

ЛЕКЦИЯ №5. ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ ФУНКЦИЙ ОДНОЙ ПЕРЕМЕННОЙ Во многих случаях применение формулы Ньютона-Лейбница для вычисления определенного интеграла затруднительно или невозможно из-за того, что первообразная оказывается слишком сложной или вообще не может быть выражена через элементарные функции. Кроме того, подынтегральная функция в задачах моделирования часто задается в табличной форме и об аналитическом представлении первообразной речь вообще не идет. В этой ситуации проблема разрешается использованием численных методов. Вначале рассмотрим интегрирование функций с использованием формул Гаусса, обладающих наивысшей алгебраической точностью. Затем приведем другие часто используемые в вычислительной практике формулы (трапеций, средних, Симпсона). 5.1. Квадратурная формула Гаусса Пусть интеграл вычисляется на стандартном интервале [1; 1] . Задача состоит в том, чтобы подобрать точки t1 ,t 2 ,...,t n и коэффициенты A1 , A2 ,..., An так, чтобы квадратурная формула 1  1 n f (t )dt   Ai f (ti ) (5.1) i 1 была точной для всех полиномов наивысшей возможной степени. Из приведенного ниже способа нахождения узлов t i и коэффициентов Ai следует, что эта наивысшая степень равна N  2n  1 . Запишем полином в виде f (t )  2 n 1 a t k 0 k k Легко показать, что для того, чтобы форму- ла (5.1) была точна для полинома, необходимо и достаточно, чтобы она была точна для степенных функций t при k  0,1,2,...,2n  1. Действительно, полагая, что согласно (5.1) k 1 n 1 i 1 k k  t dt   Ai t i k  0,1, 2,...,2n  1 (5.2) получим последовательно 1  1 1 2 n 1 2 n 1 1 2 n 1 n n 2 n 1 n 1 k 0 k 0 1 k 0 i 1 i 1 k 0 i 1 f ( t )dt    ak t k dt   ak  t k dt   ak  Ai tik   Ai  ak tik   Ai f ( ti ) , т.е. на самом деле из условия справедливости (5.2) пришли к формуле (5.1). Таким образом, система (5.2) дает 2 n соотношений для определения 2 n неизвестных Ai и t i . При этом 1 k  t dt  1 2 1  (1) k 1   , при k четном  k  1 k 1  , при k нечетном  0 Итак, согласно (5.2) коэффициенты Ai и узлы t i находятся из системы 2 n уравнений n  Ai  2,  i 1 n  Ai t i  0,  i 1 n 2 2  Ai t i  , 3  i 1 . . . . . . . . . . . . . . . . . . . . n  A t 2 n 1  0. i i  i 1   (5.3) Система (5.3) нелинейная, и ее решение найти довольно трудно. Рассмотрим следующий прием нахождения Ai и t i . Для этого нам понадобятся полиномы Лежандра, которые определяются по формуле Pn ( x)  1 dn [( x 2  1) n ] , n  0,1,2,... 2 n n! dx n Используя данную формулу, выпишем в качестве примера первые 4 полинома Лежандра P0 ( x), P1 ( x), P2 ( x), P3 ( x) P0 ( x )  1, P1 ( x )  x, P2 ( x )  1 (3 x 2  1), 2 P3 ( x )  1 (5 x 3  3 x ). 2 Полиномы Лежандра обладают рядом полезных свойств: 1) Pn (1)  1,Pn (1)  (1) n , 1 2)  P ( x) P n m n  0,1,2,... ; ( x)dx   mn N m , N m  1 2 ; 2m  1 3) полином Лежандра Pn (x) имеет n различных и действительных корней, расположенных на интервале [1;1] . (5.4) 4) Справедливо рекуррентное соотношение Pm ( x)  1 (2m  1) x Pm1 ( x)  (m  1) Pm2 ( x) . m (5.5) Теперь перейдем к рассмотрению упомянутого выше удобного способа решения нелинейной системы (5.3). Составим по узлам интегрирования многочлен n -й степени n wn ( x )   ( x  x k ) k 1 . Функция f ( x)  wn ( x) Pm ( x) при m  n  1 есть многочлен степени не выше 2n  1 . Значит для этой функции формула Гаусса справедлива: 1 n 1 i 1  wn ( x) Pm ( x) dx   Ai wn ( xi ) Pm ( xi )  0 , так как wn ( xi )  0 . Разложим wn (x ) в ряд по ортогональным многочленам Лежандра: n wn ( x)   bk Pk ( x) , k 0 1  wn ( x) Pm ( x) dx  1 1 n    bk Pk ( x)  Pm ( x)  bm N m  0 , 1 k 0 m  n 1 , bm  0 при m  n  1 . Значит wn (x ) с точностью до численного т.е. все коэффициенты множителя совпадает с Pn (x) . Таким образом, узлами формулы Гаусса являются нули многочлена Лежандра степени n . Зная t i , из линейной теперь системы первых n уравнений (5.3) легко найти коэффициенты Ai (i  1,2,..., n) . Определитель этой системы есть определитель Вандермонда. 1 Формулу  1 n f (t )dt   Ai f (t1 ) , в которой t i - нули полинома Лежандра Pn (t ) , а Ai i 1 определяют из (5.3), называют квадратурной формулой Гаусса. Пример. Вывести квадратурную формулу Гаусса для случая трех узлов, т.е. n  3 . 1. Ищем корни полинома Лежандра третьей степени P3 (t )  1 3 (5t  3t )  0 . 2 Корни полинома: t1   3 , t 2  0, t 3  5 3 . 5 2. Из первых трех уравнений (5.3) находим коэффициенты Ai A1  A2  A3  2 , 3 3 A1  A3  0 , 5 5  3 3 2 A1  A3  . 5 5 3 Отсюда 5 8 A1  A3  , A2  . 9 9 В итоге формула Гаусса при интегрировании на промежутке [1;1] имеет вид 1 1  f (t )dt  9 [5 f ( 1 3 )  8 f (0)  5 f ( 3 )] . 5 5 b При вычислении интеграла на произвольном интервале [a; b] |, т.е.  f ( x )dx для a применения квадратурной формулы Гаусса необходимо выполнить преобразование переменной x ba ba  t . 2 2 Получим f ( x)dx  ba ba ba f(  t ) dt ,  2 1 2 2 f ( x)dx  ba n  Ai f ( xi ) , 2 i 1 b  1 a тогда b  a где (5.6) xi  ba ba  ti , i  1,2,..., n ; 2 2 (5.7) здесь t i - нули полинома Лежандра Pn (t ) , т.е. Pn (t i )  0 . Погрешность формулы Гаусса с n узлами выражается формулой Rn  (b  a) 2 n1 (n!) 4 ( 2 n ) f ( ) . (2n  1)[( 2n)!]3 Отсюда, в частности, следует R3  1 b  a 7 (6) ( ) f ( ) , 15750 2 R4  1 b  a 9 (8) ( ) f ( ) 3472875 2 ......................................................... 1  b  a  (12) R8    f ( ) 648984486150  2  13 и т.д. Сделаем замечание по поводу отыскания корней полинома Лежандра произвольной степени. Эта процедура выполняется численным методом, например, можно применить метод половинного деления. Сам полином строится по рекуррентной формуле (5.5), а для начала процесса используются полиномы P0 (t ) и P1 (t ) (выписаны выше). Процедура повторяется до тех пор, пока не будут найдены все n корней полинома. При этом следует учитывать свойство полиномов (5.4), согласно которому все эти корни располагаются на интервале [1; 1] , и они все действительны и различны, т.е. кратных корней нет. 5.2. Другие формулы численного интегрирования Ставится задача вычисления определенного интеграла b F   f ( x )dx a (5.8) При построении нижеприведенных формул численного интегрирования используется общая идея, заключающаяся в том, что подынтегральную функцию f ( x ) заменяют интерполяционным многочленом, который легко интегрируется. Поскольку коэффициенты полинома линейным образом выражаются через значения интегрируемой функции в узлах, то имеет место соотношение n f ( x)   f ( xi ) i ( x)  r ( x) , (5.9) i 1 где n -количество узлов интерполяции на отрезке интегрирования [a, b] , i (x) - многочлены степени n , xi - заданные узлы интерполяции, r( x ) - остаточный член. Подставляя (5.9) в (5.8), получают формулу численного интегрирования (она называется квадратурной формулой) n F   Ai f ( xi )  R, i 0 b b a a , (5.10) Ai    i ( x)dx, R   r ( x)dx где Ai -веса квадратурной формулы, а R - погрешность или остаточный член формулы. Понятно, что формула (5.10) точна для полинома степени n , а следовательно и для всех степеней x от 0 до n , т.е. справедливы соотношения, использованные ранее при получении системы (5.3) b n a i 0 k k  x dx   Ai xi при k  0,1, 2....n . При этом b k  x dx  a b k 1  a k 1 . k 1 (5.11) Отсюда, задаваясь количеством и расположением узлов можно, применяя (5.11), вычислить веса квадратурной формулы. Например, при n  0, x0  ab получим b  a  A0 , т.е. квадратурная формула име2 ет вид b  f ( x)dx  (b  a) a f( ab ) 2 (5.12) Эта формула называется формулой средних. При n  1, x0  a, x1  b получим формулу трапеций. При n  2, x 0  a, x1  b  f ( x)dx  a ab , x 2  b будет построена формула Симпсона 2 (b  a)  f ( x0 )  4 f ( x1 )  f ( x2 ) . 6 (5.13) Эту процедуру можно продолжить, получая новые формулы со все большим количеством узлов. Однако можно поступить иначе. Используя интерполяционный полином Лагранжа, т.е. подставляя в (5.9) в качестве функций  i ( x ) лагранжевы коэффициенты L(i n ) ( x ) , и проведя интегрирование в формуле (5.10), вычислим Ai Ai  ( b  a )H i , где H i - коэффициенты Котеса. Существуют таблицы коэффициентов Котеса вплоть до значений n=10 Приведем сводку ряда простейших формул, используемых в практике численного интегрирования, при большом количестве узлов и постоянном расстоянии между узлами (шаге сетки). Формула трапеций: b 1  f ( x)dx  h ( 2 f a  f1 ... f N 1  1 fN ) , 2 b 1 R   h 2  f ( x) dx  O(h 2 ) , 12 a h  xi  xi1  const , где R - асимптотическая погрешность формулы; f i  f ( xi ) - значения интегрируемой функции в узлах. Формула средних: b N  f ( x)dx  h f ( x i 1 a i 1 / 2 ), xi 1 / 2  xi 1  xi , 2 b 1 2 R h  f ( x)dx  O(h 2 ) . 24 a Формула Симпсона: b h N 1 2  f ( x)dx  3  ( f i 0 a 2i  4 f 2i 1  f 2i 2 ) , b 1 4 ( 4) R h  f ( x)dx  O(h 4 ) . 180 a Отметим, что если подинтегральная функция не имеет соответствующих производных, то указанный теоретический порядок точности не достигается. Так, если на отрезке интегрирования не существуют 3-я и 4-я производные, то порядок точности формулы Симпсона будет только 2-ой, O(h2).