1611688890-f641c9ec8276824e4686da772eb56520 (826652), страница 27
Текст из файла (страница 27)
2.23), коль скоро основано на приближении подинтегральной функции подходящей параболой.20Читатель может самостоятельно убедиться, что та же самая формула получается (но только более длинно и утомительно) в результатеинтегрирования по [a, b] интерполяционного полинома второй степенив форме Лагранжаa + ba + b(x − b)x−(x − a)(x − b)2f (a) + fP2 (x) = a+ba+ba+b2(a − b)−a−ba−222a + b(x − a) x −2+f (b)a + b(b − a) b −2a + b2a + b=(x−b)f(a)−2(x−a)(x−b)fx−(b − a)222a + bf (b) ,+ (x − a) x −220 Приведённый нами элегантный вывод формулы Симпсона восходит к учебникуА.Н. Крылова [18].1602.
Численные методы анализакоторый строится для подинтегральной функции по узлам a, (a + b)/2и b.Предложение 2.12.1 Квадратурная формула Симпсона имеет алгебраическую степень точности 3, т. е. является точной для любогополинома степени не выше третьей.Доказательство. Отметим прежде всего, что для полиномов степенине выше 2 точность квадратурной формулы Симпсона следует прямо изтого, что она построена как интерполяционная квадратурная формула, в которой подинтегральная функция интерполируется полиномомвторой степени.
Поэтому достаточно показать, что формула Симпсонаточна для монома x3 , но не является точной для более высоких степеней.При интегрировании x3 получаемZbx3 dx =ab 4 − a4.4С другой стороны, согласно формуле Симпсона a + b 3b−a 3b − a 3 a3 + 3a2 b + 3ab2 + b333a +4a ++b =+b6262==b − a 3a3 + 3a2 b + 3ab2 + 3b3·62 b 4 − a4b−a 3a + a2 b + ab2 + b3 =,44что совпадает с результатом точного интегрирования.Для монома x4 длинными, но несложными выкладками нетруднопроверить, что результат, даваемый формулой Симпсона для интеграла по интервалу [a, b], т.
е. a + b 4b−a 4a +4+ b4 ,62отличается от точного значения интегралаZabx4 dx =b 5 − a551612.12. Численное интегрированиена величину (b − a)5 /120. Она не зануляется при a 6= b, так что наполиномах четвёртой степени формула Симпсона уже не точна.Итак, несмотря на то, что формула Симпсона основана на интерполяции подинтегральной функции полиномом степени 2, фактическаяточность формулы более высока, чем та, что обеспечивается полиномомвторой степени. В этой ситуации, желая аккуратно оценить погрешность формулы Симпсона на основе известной погрешности алгебраической интерполяции (аналогично выводу погрешности формулы трапеций в §2.12б), желательно в выражении для погрешности взять болеевысокую степень переменной.
Иными словами, при оценке погрешности формулы Симпсона нужно воспользоваться для подинтегральнойфункции интерполяционным полиномом третьей степени, а не второй,на основе которого она была построена. При наличии всего трёх узлов мы находимся тогда в условиях задачи интерполяции с кратнымиузлами.Предполагая существование производной f ′ в среднем узле x1 =(a + b)/2, можно считать, к примеру, что именно он является кратнымузлом. Формально нам необходим такой интерполяционный полином3-й степени H3 (x), чтоH3 (a) = f (a),H3a + b2=fa + b2(2.121)H3 (b) = f (b),,H3′a + b2= f′a + b2,(2.122)хотя конкретное значение производной в средней точке (a + b)/2 далееникак не будет использоваться.
Здесь нам важно лишь то, что при любом значении этой производной решение задачи интерполяции (2.121)–(2.122) существует, и потребуется оценка его погрешности.Существование и единственность решения подобных задач былаустановлена в §2.4 и там же обосновывается оценка его погрешности(2.45): nf (m+1) ξ(x) Y(2.123)f (x) − Hm (x) =(x − xi )Ni ,(m + 1)! i=0где Ni — кратности узлов, m = N0 + N1 + . . . + Nn − 1 — степень интерполяционного полинома, а ξ(x) — некоторая точка из [a, b], зависящая1622. Численные методы анализаот x. Для решения задачи (2.121)–(2.122) справедливоf (4) ξ(x)a + b 2f (x) − H3 (x) =(x − b).· (x − a) x −242Далее, из того, что формула Симпсона точна для полиномов третьей степени, а также из условий (2.121)–(2.122) следуют равенстваZ ba + bb−aH3 (x) dx =+ H3 (b)· H3 (a) + 4H362aa + bb−a+ f (b) .(2.124)· f (a) + 4f=62Отсюда уже нетрудно вывести выражение для погрешности квадратурной формулы Симпсона:Z ba + bb−af (x) dx −R(f ) =+ f (b)· f (a) + 4f62a=ZbZba=af (x) − H3 (x) dxв силу (2.124)f (4) ξ(x)a + b 2· (x − a) x −(x − b) dx242из (2.123).Из него следует оценкаZ b f (4) ξ(x)a + b 2|R(f )| = (x − b) dx · (x − a) x − a242≤ZbaM4≤24 f (4) ξ(x) a + b 2 (x − b) dx · (x − a) x −242Z ba + b 2(x − b) dx ,· (x − a) x − a2поскольку в интегрируемой функции подвыражениеa + b 2(x − a) x −(x − b)2(2.125)1632.12.
Численное интегрированиене меняет знак на интервале интегрирования [a, b]. В (2.125), как обычно, обозначено M4 = maxx∈[a,b] |f (4) (x)|.Для вычисления фигурирующего в (2.125) интеграла сделаем замену переменныхa+bt=x−,2тогдаZ ba + b 2(x − b) dx(x − a) x −2a==Zb−a2Zb−a2b − a 2b − at+t t−dtb−a22−2b−a− 2(b − a)2(b − a)5t2 t2 −dt = −.4120ОкончательноM4 (b − a)5.2880Как видим, более тонкие рассуждения о свойствах формулы Симпсонапозволили получить действительно более точную оценку её погрешности.|R(f )| ≤2.12гОбщие интерполяционныеквадратурные формулыКвадратурными формулами интерполяционного типа мы назвали(см. §2.12б) формулы, получающиеся в результате замены подинтегральной функции f (x) интерполяционным полиномом Pn (x), которыйпостроен по некоторой совокупности простых узлов x0 , x1 , .
. . , xn изинтервала интегрирования. Выпишем для общего случая этот полиномв форме Лагранжа:Pn (x) =nXf (xi ) φi (x),i=0гдеφi (x) =(x − x0 ) · · · (x − xi−1 )(x − xi+1 ) · · · (x − xn+1 ),(xi − x0 ) · · · (xi − xi−1 )(xi − xi+1 ) · · · (xi − xn+1 )1642. Численные методы анализа— базисные полиномы Лагранжа (стр. 54).Интерполяционная квадратурная формула должна получаться изприближёного равенстваZ bZ bf (x) dx ≈Pn (x) dx(2.126)aaв результате выполнения интегрирования в правой части. Как следствие, в представлении (2.115) весовые коэфициенты формулы имеютвидZ bci =φi (x) dxa=Z(2.127)ba(x − x0 ) · · · (x − xi−1 )(x − xi+1 ) · · · (x − xn+1 )dx,(xi − x0 ) · · · (xi − xi−1 )(xi − xi+1 ) · · · (xi − xn+1 )i = 0, 1, .
. . , n. Эти значения весов ci , однозначно определяемые по узлам x0 , x1 , . . . , xn , являются отличительным характеристическим признаком именно интерполяционной квадратурной формулы. Если длязаданного набора узлов у какой-либо квадратурной формулы весовыекоэффициенты равны (2.127), то можно считать, что она построенана основе алгебраической интерполяции подинтегральной функции поэтим узлам, взятым с единичной кратностью.Теорема 2.12.1 Для того, чтобы квадратурная формула (2.115), построенная по (n + 1) попарно различным узлам, была интерполяционной, необходимо и достаточно, чтобы её алгебраическая степеньточности была не меньшей n.В качестве замечания к формулировке нужно отметить, что в условиях теоремы квадратурная формула на самом деле может иметь алгебраическую степень точности выше n, как, например, формула средних прямоугольников или формула Симпсона.Доказательство.
Необходимость условий теоремы очевидна: интерполяционная квадратурная формула на n + 1 узлах, конечно же, точнана полиномах степени n, поскольку тогда подинтегральная функциясовпадает со своим алгебраическим интерполянтом.Покажем достаточность: если квадратурная формула (2.115), построенная по (n + 1) узлу, является точной для любого алгебраического полинома степени n, то её весовые коэффициенты вычисляются по1652.12. Численное интегрированиеформулам (2.127), т. е. она является квадратурной формулой интерполяционного типа.В самом деле, для базисных интерполяционных полиномов φi (x)выполнено свойство (2.9)(0, при i 6= j,φi (xj ) = δij =1, при i = j,и они имеют степень n. Следовательно, применяя рассматриваемуюквадратурную формулу для вычисления интеграла от φi (x), получимZbφi (x) dx =anXck φi (xk ) =k=0nXck δik = ci ,k=0и это верно для всех i = 0, 1, .
. . , n. Иными словами, имеет место равенство (2.127), что и требовалось доказать.В частности, если в интерполяционной квадратурной формуле вместо подинтегральной функции взять полином P0 (x) = x0 = 1, то получаем равенствоZ bnX1 dx =ck ,b−a =ak=0— сумма весов такой квадратурной формулы равна длине интервалаинтегрирования.Из (2.126) ясно, что погрешность интерполяционных квадратурныхформул равнаZ bRn (f, x) dx,R(f ) =aгде Rn (f, x) — остаточный член алгебраической интерполяции.
В §2.2дбыла получена оценка для Rn (f, x) в форме Коши (2.25)f (n+1) ξ(x)· ωn (x),Rn (f, x) =(n + 1)!где ξ(x) ∈ [a, b], и поэтому1R(f ) =(n + 1)!Zbaf (n+1) ξ(x) ωn (x) dx.1662. Численные методы анализаСправедлива огрублённая оценка|R(f )| ≤Mn+1(n + 1)!Zab(2.128)|ωn (x)| dx,где Mn+1 = maxx∈[a,b] |f (n+1) (x)|. Из неё можно ещё раз заключить,что квадратурная формула интерполяционного типа, построенная по(n + 1) узлам, является точной для любого полинома степени не болееn, поскольку тогда Mn+1 = 0.Из наших рассуждений видно, что оценка (2.128) является простейшей, использующей лишь основные свойства алгебраического интерполянта. В некоторых случаях она может оказаться существенно завышенной, что мы могли видеть на примере формулы Симпсона.2.12дДальнейшие формулы Ньютона-КотесаВ §2.12б и §2.12в простейшие квадратурные формулы НьютонаКотеса — формулы прямоугольников и трапеций, формула Симпсона— были выведены и исследованы средствами, индивидуальными длякаждой отдельной формулы.
В этом разделе мы взглянем на формулыНьютона-Котеса с более общих позиций.Зафиксировав номер n, n ≥ 1, возьмём на интервале интегрирования [a, b] равноотстоящие друг от друга узлы(n)xk= a + kh,k = 0, 1, . . . , n,h=b−a.nДля определения весов формул Ньютона-Котеса необходимо вычислить величины (2.127), которые мы обозначим для рассматриваемогочастного случая как(n)Ak=Zab(n) (n) (n) · · · x − xk−1 x − xk+1 · · · x − xn+1dx,(n)(n) (n) (n)(n) (n)(n) (n)xk − x0 · · · xk − xk−1 xk − xk+1 · · · xk − xn(n) x − x0k = 0, 1, . .
. , n. Сделаем в этом интеграле замену переменных x = a+th,1672.12. Численное интегрированиегде t пробегает интервал [0, n]. Тогдаdx = h dt,(n) (n) (n) (n) · · · x − xk−1 x − xk+1 · · · x − xn+1x − x0= hn t(t − 1) · · · (t − k + 1)(t − k − 1) · · · (t − n),(n)(n) (n)(n)(n) (n) (n)(n) xk − x0 · · · xk − xk−1 xk − xk+1 · · · xk − xn= (−1)n−k hn k!(n − k)!,где считается, что 0! = 1. Окончательно(n)Ak(−1)n−k=hk!(n − k)!Znt(t − 1) · · · (t − k + 1)(t − k − 1) · · · (t − n) dt,0k = 0, 1, . . . , n. Чтобы придать результату не зависящий от интервалаинтегрирования вид, положим(n)Ak(n)= (b − a) Bk ,где(n)Bk(−1)n−k=n k!(n − k)!Znt(t − 1) · · · (t − k + 1)(t − k − 1) · · · (t − n) dt.0(n)Теперь уже величины Bkкоэффициентами Котеса.К примеру, для n = 1Z(1)B0 = −(1)B1=10Z01не зависят от h и [a, b].
Они называются(t − 1) dt = −11(t − 1)2 = 2,201t2 1t dt == .2 02Мы вновь получили веса квадратурной формулы трапеций (2.117). Для1682. Численные методы анализаслучая n = 2Z(2)B01=4(2)B11=−2(2)B1 =1420Z1(t − 1)(t − 2) dt =40Z021t(t − 2) dt = −22t(t − 1) dt =1421t3t2− 3 + 2t = ,326024t32 −t = ,36021t3t2 = .−32 06Полученные коэффициенты соответствуют формуле Симпсона (2.120).И так далее.За прошедшие три с лишним столетия коэффициенты Котеса былитщательно вычислены для значений n из начального отрезка натурального ряда. В Табл. 2.2, заимствованной из книги [18], приведены коэффициенты Котеса для n ≤ 10 (см.