1611688890-f641c9ec8276824e4686da772eb56520 (826652), страница 26
Текст из файла (страница 26)
То, насколько широким являетсякласс функций, на котором точна рассматриваемая формула, можетслужить косвенным признаком её точности вообще. Очень часто в качестве класса «пробных функций» F, для которых исследуется совпадение результата квадратурной формулы и искомого интеграла, беруталгебраические полиномы. В этой связи полезноОпределение 2.12.1 Алгебраической степенью точности квадратурной формулы называют наибольшую степень алгебраических полиномов, для которых эта квадратурная формула является точной.Сответственно, из двух квадратурных формул более предпочтительной можно считать ту, которая имеет бо́льшую алгебраическуюстепень точности.
Неформальным обоснованием этого критерия служит тот факт, что с помощью полиномов более высокой степени можно получать более точные приближения функций, как локально (с помощью формулы Тейлора), так и глобально (к примеру, с помощьюразложения по полиномам Чебышёва или Лежандра).Рассмотрим теперь влияние погрешностей реальных вычислений наответ, получаемый с помощью квадратурных формул.
Предположим,что значения f (xk ) интегрируемой функции в узлах xk вычисляютсянеточно, с погрешностями δk . Тогда при вычислениях по квадратурнойформуле получимnXk=0ck f (xk ) + δk=nXk=0ck f (xk ) +nXck δk .k=0Если для всех k = 0, 1, . . . , n знаки погрешностей δk совпадают со знаками весов ck , то общая абсолютная погрешность результата,полученPного по квадратурной формуле, становится равной k |ck | |δk |.
Желая1522. Численные методы анализанайти оценку сверху для этого выражения, можем вынести maxk |δk |,как общий множитель, за знак суммы. Следовательно, величинуnXk=0|ck |,т. е. сумму модулей весов квадратурной формулы, имеет смысл рассматривать как коэффициент усиления погрешности при вычисленияхс этой формулой.Если при значительном количестве узлов n мы хотим организоватьвычисления по квадратурной формуле наиболее устойчивым образом,то все весовые коэффициенты ck должны быть положительны: именнотогда при прочих равных условиях сумма модулей весов минимальна.В частности, в случае интегрирования функций, принимающих значения одного знака, мы избегаем тогда потери точности при вычитанииблизких значений, которое могло бы случиться в формуле, где одновременно присутствуют положительные и отрицательные веса.2.12бФормулы Ньютона-Котеса.Простейшие квадратурные формулыПростейший приём построения квадратурных формул — замена подинтегральной функции f (x) на интервале интегрирования [a, b] на «более простую», легче интегрируемую функцию, которая интерполируетили приближает f (x) по заданным узлам x0 , x1 , .
. . , xn . В случае,когда f (x) заменяется интерполянтом и все рассматриваемые узлы —простые, говорят о квадратурных формулах интерполяционного типа,или, что равносильно, об интерполяционных квадратурных формулах.Наиболее часто подинтегральную функцию интерполируют полиномами, и в нашем курсе мы будем рассматривать только такие интерполяционные квадратурные формулы.Формулами Ньютона-Котеса называют интерполяционные квадратурные формулы, полученные с помощью алгебраической интерполяции подинтегральной функции на равномерной сетке с простыми узлами. В зависимости от того, включаются ли концы интервала интегрирования [a, b] в множество узлов квадратурной формулы или нет,различают формулы Ньютона-Котеса замкнутого типа и открытоготипа.2.12.
Численное интегрирование153Рис. 2.21. Иллюстрация квадратурных формуллевых и правых прямоугольниковДалее мы построим и исследуем формулы Ньютона-Котеса для n =0, 1, 2, причём будем строить наиболее популярные формулы замкнутого типа за исключением случая n = 0, когда имеем всего один узели замкнутая формула просто невозможна.Если n = 0, то подинтегральная функция f (x) интерполируется полиномом нулевой степени, т. е. какой-то константой, равной значениюf (x) в единственном узле x0 .
Если взять x0 = a, то получается квадратурная формула «левых прямоугольников», а если x0 = b — формула«правых прямоугольников» (см. Рис. 2.21).Ещё один естественный вариант выбора единственного узла —x0 = 12 (a + b),т. е. как середины интервала интегрирования [a, b]. При этом приходимк квадратурной формулеZ ba + b,f (x) dx ≈ (b − a) · f2aназываемой формулой средних прямоугольников: согласно ей интегралберётся равным площади прямоугольника с основанием (b − a) и высотой f (a + b)/2 (см. Рис. 2.22). Эту формулу нередко называют также просто «формулой прямоугольников», так как она является наиболее часто используемым вариантом рассмотренных простейших квадратурных формул.Оценим погрешность формулы средних прямоугольников методомлокальных разложений, который ранее был использован при исследовании численного дифференцирования.
Разлагая f (x) в окрестности1542. Численные методы анализаточки x0 = 12 (a + b) по формуле Тейлора с точностью до членов первого порядка, получимf (x) = fa + b2+f ′a + b a + b f ′′ (ξ) a + b 2· x−+, (2.116)· x−2222где ξ — зависящая от x точка интервала [a, b], которую корректно обозначить через ξ(x). Остаточный член равенR(f ) =Zbf (x) dx − (b − a) · faa + b2Z b a + b dx=f (x) − f2a=Z b a + b a + b f ′′ (ξ(x)) a + b 2f′· x−+dx· x−2222a=Zaпосколькуbf ′′ (ξ(x)) a + b 2dx,· x−22Z baa + bx−dx =2Zb−a2b−a− 2t dt = 0,— интеграл от первого члена разложения (2.116) зануляется.
Следовательно, с учётом принятого нами ранее обозначенияMp := max f (p) (x)x∈[a,b]можно выписать оценкуZ bZ b ′′ 2 f (ξ) a + b 2 · x − a + b dx ≤ M2dxx−|R(f )| ≤ 2 22 a2aba + b 3 M2 1 M2 (b − a)3x−==·.2 3224aОтсюда, в частности, следует, что для полиномов степени не выше 1формула (средних) прямоугольников даёт точное значение интеграла,1552.12. Численное интегрированиеРис. 2.22.
Иллюстрация квадратурных формулсредних прямоугольников и трапецийколь скоро вторая производная подинтегральной функции тогда зануляется и M2 = 0.Полученная оценка точности неулучшаема, так как достигается на2функции g(x) = x − 21 (a + b) . При этомM2 = max |g ′′ (x)| = 2,gx∈[a,b]a + b2= 0,и потомуZabg(x) dx − (b − a) · ga + b2=(b − a)3M2 (b − a)3=,1224т. е. имеем точное равенство на погрешность.Рассмотрим теперь квадратурную формулу Ньютона-Котеса, соответствующую случаю n = 1, когда подинтегральная функция приближается интерполяционным полиномом первой степени.
Для формулызамкнутого типа построим его по узлам x0 = a и x1 = b, совпадающимс концами интервала интегрирования:P1 (x) =x−ax−bf (a) +f (b).a−bb−a1562. Численные методы анализаИнтегрируя это равенство, получимZ bZZf (a) bf (b) bP1 (x) dx =(x − b) dx +(x − a) dxa−b ab−a aabbf (a) (x − b)2 f (b) (x − a)2 = + b−aa−b22aab−af (a) + f (b) .2Мы вывели квадратурную формулу трапецийZ bf (a) + f (b)f (x) dx ≈ (b − a) ·,2a=(2.117)название которой также навеяно геометрическим образом. Фактически, согласно этой формуле точное значение интеграла заменяется назначение площади трапеции (стоящей боком на оси абсцисс) с высотой(b − a) и основаниями, равными f (a) и f (b) (см. Рис.
2.22).Чтобы найти погрешность формулы трапеций, вспомним оценку(2.25) для погрешности интерполяционного полинома. Из неё следует,чтоf ′′ (ξ(x))f (x) − P1 (x) =· (x − a)(x − b)2для некоторой точки ξ(x) ∈ [a, b]. Таким образом, для формулы трапеций остаточный член естьZ b ′′Z bf (ξ(x))· (x − a)(x − b) dx,f (x) − P1 (x) dx =R(f ) =2aaно вычисление полученного интеграла на практике нереально из-занеизвестного вида ξ(x). Как обычно, имеет смысл вывести какие-то более удобные оценки погрешности, хотя они, возможно, будут не стольточны.Поскольку выражение (x − a)(x − b) всюду на интервале [a, b] кромеего концов сохраняет один и тот же знак, тоZ b ′′f (ξ(x))|R(f )| ≤· |(x − a)(x − b)| dx2aZM2 b≤· (x − a)(x − b) dx ,2a1572.12. Численное интегрированиегде M2 = maxx∈[a,b] |f ′′ (x)|.
ДалееZab(x − a)(x − b) dx =Zbax2 − (a + b) x + ab dxbbbx3 x2 =− (a + b) + abx 3 a2 aa==(2.118)1 32 b − a3 − 3(a + b) b2 − a2 + 6ab(b − a)61(b − a)3−b3 + 3ab2 − 3a2 b + a3 = −.66Поэтому окончательно|R(f )| ≤M2 (b − a)3.12Эта оценка погрешности квадратурной формулы трапеций неулучшаема, поскольку достигается при интегрировании функции g(x) = (x−a)2по интервалу [a, b].2.12вКвадратурная формула СимпсонаПостроим квадратурную формулу Ньютона-Котеса для n = 2, т. е.для трёх равномерно расположенных узловx0 = a,x1 =a+b,2x2 = bиз интервала интегрирования [a, b].Для упрощения рассуждений выполним параллельный перенос криволинейной трапеции, площадь которой мы находим с помощью интегрирования, и сделаем точку a началом координат оси абсцисс (см.Рис.
2.23). Тогда правым концом интервала интегрирования станетl = b − a. ПустьP̆2 (x) = c0 + c1 x + c2 x2— полином второй степени, интерполирующий сдвинутую подинтегральную функцию по узлам 0, l/2 и l. Если график P̆2 (x) проходит1582. Численные методы анализа99Kab0lРис. 2.23. Иллюстрация вывода квадратурной формулы Симпсоначерез точки плоскости с координатами0, f (a) ,тоl a + b,,f22l, f (b) ,c0 = f (a),a + bll2c0 + c1 + c2 = f,242c0 + c1 l + c2 l2 = f (b).(2.119)Площадь, ограниченная графиком интерполяционного полинома, равнаZ0ll3l2c0 + c1 x + c2 x2 dx = c0 l + c1 + c223=l6c0 + 3c1 l + 2c2 l2 .6Фактически, для построения квадратурной формулы требуется решитьотносительно c0 , c1 и c2 систему уравнений (2.119) и потом подставитьрезультаты в полученное выше выражение.
Но можно выразить трёхчлен 6c0 + 3c1 l + 2c2 l2 через значения подинтегральной функции f вузлах, не решая систему (2.119) явно.2.12. Численное интегрирование159Умножая второе уравнение системы (2.119) на 4 и складывая с первым и третьим уравнением, получимa + b6c0 + 3c1 l + 2c2 l2 = f (a) + 4f+ f (b).2Таким образом,Z lZ lc0 + c1 x + c2 x2 dxP̆2 (x) dx =00=l6a + bf (a) + 4f+ f (b) ,2что даёт приближённое равенствоZ ba + bb−af (x) dx ≈+ f (b) .· f (a) + 4f62a(2.120)Оно называется квадратурной формулой Симпсона или формулой парабол (см. Рис.