Блюмин А.Г., Федотов А.А., Храпов П.В. - Численные методы (1078410), страница 2
Текст из файла (страница 2)
Пусть f ∈ C 1 [a, b], ξi ∈ [a, b] — произвольные точки, i = 1, 2, ..., n.Тогда существует такая тока ξ ∈ [a, b], что[f (ξ1 ) + f (ξ2 ) + ... + f (ξn )]/n = f (ξ).Эта лемма вытекает из очевидных неравенствmin f (x) 6 [f (ξ1 ) + f (ξ2 ) + ... + f (ξn )]/n 6 max f (x)x∈[a,b]x∈[a,b]и теоремы о промежуточных значениях непрерывной функции.Используя лемму, получаем формулу средних прямоугольников с остаточным членомhZ2f (x)dx = h · f0 +h3 00h· f (ξ), | ξ |6 .242− h29Глава 1. Численные методы вычисления определенного интеграла1.2.3. Формула трапецийПусть f ∈ C 2 [0, h].
ПолагаемZhf (x)dx ≈ h ·I=f0 + f1,2(1.2.4)0где f0 = f (0), f1 = f (h). Из формулы (1.2.4) видно, что искомое значение интеграла приближенно заменяется величиной площади закрашенной нарис. (1.1,б ) трапеции.Аналогично тому, как это сделано в п. (1.2.2) можно получить формулутрапеций с остаточным членомZhf0 + f1 h3 00−· f (ξ), ξ ∈ [0, h].f (x)dx = h ·21201.2.4. Формула СимпсонаПредположим, что f ∈ C 4 [−h, h] и требуется вычислить интегралZhI=f (x)dx.−hЗначение этого интегралаyприближенно заменяем величи-p(x)f(x)ной площади закрашенной криволинейной трапеции (рис. 1.2),f1ограниченной сверху параболойf0f-1p(x), проходящей через точки(−h, f−1 ), (0, f0 ),(h, f1 ), гдеfi = f (i · h), i = −1, 0, 1.
Эта-hпарабола задается уравнением0Рис. 1.2.10hxГлава 1. Численные методы вычисления определенного интегралаp(x) = f0 +иf1 − f−1f−1 − 2f0 + f1 2·x+·x2h2h2Zhp(x)dx =h· (f−1 + 4f0 + f1 ).3f (x)dx ≈h· (f−1 + 4f0 + f1 ).3−hСледовательноZh(1.2.5)−hФормула Симпсона с остаточным членом имеет видZhh5 (IV )h·f(ξ), ξ ∈ [−h, h].f (x)dx = · (f−1 + 4f0 + f1 ) −390−hРассмотренные квадратурные формулы средних прямоугольников (1.2.2),трапеций (1.2.4) и Симпсона (1.2.5) назовем каноническими.1.2.5. Составные квадратурные формулыНа практике, если требуется вычислить приближенно интеграл, обычноделят заданный отрезок [a, b] на n равных частичных отрезков [xi−1 , xi ], гдеxi = a + i · h, i = 0, 1, ..., n; x0 = a, xn = b, h = (b − a)/n.
На каждом частичном отрезке используют каноническую квадратурную формулу и суммируютполученные результаты. При применении формул средних прямоугольникови трапеций длину частичных отрезков удобно принять за h, а при использовании формулы Симпсона — за 2h. В результате получаются следующиеформулы, которые будем называть составными.11Глава 1. Численные методы вычисления определенного интегралаСоставная квадратурная формула средних прямоугольников записывается в виде (рис. 1.3)Zbf (x)dx ≈ h · (fc1 + fc2 + ...
+ fcn ),(1.2.6)aгде h = (b − a)/n, fci = f (xci ); xci = a + (i − 1/2)h, i = 1, 2, ..., n — координатысредних точек частичных отрезков [xi−1 , xi ].yf(x)f(b)f(a)a xc1xc2...xci...xcn bxРис. 1.3.Погрешность Rn получается в результате суммирования погрешностей почастичным отрезкамnh3 X 00h3Rn =f (ξi ) = n ·24 i=124!n1 X 00f (ξi ) ,n i=1где xi−1 < ξi < xi . В соответствии со сформулированной выше леммой последнее выражение для Rn можно переписать в видеh3(b − a) 00Rn =· n · f 00 (ξ) = h2 ·· f (ξ), ξ ∈ [a, b].242412Глава 1. Численные методы вычисления определенного интегралаПусть M — максимальное значение модуля второй производной функцииf (x) на отрезке [a, b], т.е. M = max | f 00 (x) |; тогда из выражения для Rnx∈[a,b]получаем следующую оценку:| Rn |6 h2 ·(b − a) · M,24это означает, что погрешность формулы средних прямоугольников на всемотрезке интегрирования [a, b] есть величина O(h2 ) (см.
определение 2).В этом случае говорят, что квадратурная формула имеет второй порядокточности.Замечание. Возможны формулы прямоугольников и при ином, чем в формуле средних прямоугольников, расположении узлов. Например,Zbf (x)dx ≈nXZbh · fxi−1 ,f (x)dx ≈i=1aanXh · fxi .i=1Однако из-за нарушения симметрии погрешность таких формул является величиной O(h), т.е. порядок точности таких формул на единицу ниже порядкаточности формулы средних прямоугольников.Составная квадратурная формула трапеций имеет видZbf (x)dx ≈ h ·fnf0+ f1 + f2 + ... + fn−1 +22,(1.2.7)aгде fi = f (xi ), xi = a + i · h, h = (b − a)/n, i = 0, 1, ..., n.Аналогично предыдущему случаю можно получить выражение для погрешности Rn составной формулы трапецийRn = −h2 ·(b − a) 00· f (ξ), ξ ∈ [a, b].12Тогда имеет место оценка| Rn |6 h2 ·(b − a) · M, M = max | f 00 (x) | .12x∈[a,b]13Глава 1.
Численные методы вычисления определенного интегралаТаким образом, формула трапеций (1.2.7) имеет, так же как и формула средних прямоугольников (1.2.6), второй порядок точности (Rn = O(h2 ));следует заметить, что ее погрешность оценивается величиной в два раза большей, чем погрешность формулы средних прямоугольников.Составная квадратурная формула Симпсона записывается такZbf (x)dx ≈h·3f0 + f2n + 4 ·nXf2i−1 + 2 ·i=1an−1X!f2i ,(1.2.8)i=1где fj = f (xj ), xj = a + j · h, h = (b − a)/(2n), j = 0, 1, ..., 2n.Погрешность составной формулы Симпсона имеет видRn = −h4 ·(b − a) (IV )·f(ξ), ξ ∈ [a, b].180Отсюда получаем оценку| Rn |6 h4 ·(b − a) · M, M = max | f (IV ) (x) |,180x∈[a,b]т.е. составная формула Симпсона существенно точнее, чем формулы среднихпрямоугольников и трапеций.
Она имеет на отрезке [a, b] четвертый порядокточности (Rn = O(h4 )).Из выражений погрешностей видно, что формулы средних прямоугольников и трапеций точны для многочленов первой степени, т.е. для линейныхфункций, а формула Симпсона точна для многочленов третьей степени (дляних погрешность равна нулю).1.2.6. Квадратурные формулы ГауссаБудем считать, что интеграл предварительно приведен к стандартнойформе, когда областью интегрирования является отрезок [−1, 1]. Итак, пусть14Глава 1. Численные методы вычисления определенного интегралатребуется вычислить интегралZ1I=f (x)dx.−1Мы рассматривали до сих пор квадратурные формулы с заданными узлами и убедились, что формулы средних прямоугольников и трапеций точныдля многочленов первой степени, а формула Симпсона — для многочленовтретьей степени.
Пусть мы имеем квадратурную формулу с n узловыми точкамиIn =nXqi · f (xi ).(1.2.9)i=1Если считать неизвестными не только весовые коэффициенты qi , но и узлыxi , то можно потребовать, чтобы квадратурная формула (1.2.9) была точна для полиномов наиболее высокой степени m.
Такую формулу называютквадратурной формулой Гаусса. При этом оказывается, что m = 2n − 1.Формула (1.2.9) должна быть точна для f (x) = 1, x, x2 , ..., x2n−1 , т.е.In =nXi=1qi · xli =Z1−1ll+1 1x = 1 + (−1) ,xl dx =l + 1 −1l+1где l = 0, 1, ..., 2n−1. В результате для узлов xi и коэффициентов qi получимследующую систему 2n нелинейных уравнений:q1 + q2 + ... + qn = 2 ,q1 x1 + q2 x2 + ...
+ qn xn = 0 ,2q1 x21 + q2 x22 + ... + qn x2n = ,3...2n−1 q1 x2n−1 + q2 x2n−1 + ... + qn x2n−1 = 1 + (−1).n122n15(1.2.10)Глава 1. Численные методы вычисления определенного интегралаВ простейшем случае n = 1 систему (1.2.10) можно решить и убедиться,в том, что полученная формула Гаусса совпадает с формулой средних прямоугольников: In = 2 · f (0), и что она верна для любой линейной функцииf (x) = c0 + c1 x. В общем случае при произвольном n можно показать (см.,например, [4]), что узлами квадратурной формулы Гаусса являются корниполинома Лежандра Pn (x), а весовые коэффициенты вычисляются по формулеZ1qj =Qn−1,j (x)dx, j = 1, ..., n;(1.2.11)−1где подынтегральная функцияQn−1,j (x) =(x − x1 )(x − x2 ) · · · (x − xj−1 )(x − xj+1 ) · · · (x − xn ).(xj − x1 )(xj − x2 ) · · · (xj − xj−1 )(xj − xj+1 ) · · · (xj − xn )Функция Qn−1,j (x) является полиномом степени (n − 1).
В числителе унего стоит произведение (n − 1)-ого множителей (x − xi ), i = 1, ..., n, i 6= j;в знаменателе — значение числителя в узле x = xj . Таким образом, полиномQn−1,j (x) в узлах xi принимает следующие значения 0, i 6= j,Qn−1,j (xi ) = 1, i = j.Полиномы Лежандра. Они определяются формулойPn (x) =1 dn 2(x − 1)n , n = 0, 1, 2, ...,nn2 n! dx(1.2.12)Согласно (1.2.12) P0 (x)=1, P1 (x) = x. Для последующих значений n можновоспользоваться рекуррентным соотношениемnPn (x) = (2n − 1) · x · Pn−1 (x) − (n − 1)Pn−2 (x).16Глава 1. Численные методы вычисления определенного интегралаПользуясь этой формулой, выпишем полиномы Лежандра для n = 2, 3, 4, 5:1P2 (x) = (3x2 − 1) ,21P4 (x) = (35x4 − 30x2 + 3) ,81P3 (x) = (5x3 − 3x) ,21P5 (x) = (63x5 − 70x3 + 15x).8Графики полиномов Pn (x) до n = 5 представлены на рис. 1.4.aPn(x)Pn(x)P2(x)P5(x) P3(x) 0.500.50.501.0 x-0.5ɛ0.5-0.5-0.5P0(x)1.0xP4(x)P1(x)-1.0-1.0Рис. 1.4.Полиномы Лежандра с четными номерами являются четными функциями, а с нечетными — нечетными.
Полиномы Лежандра Pn (x) в точках x = ±1принимают следующие значения: Pn (1) = 1, Pn (−1) = (−1)n . На интервале(−1, 1) многочлен Pn (x) имеет n простых нулей. В силу четности или нечетности Pn (x) нули полиномов Лежандра располагаются симметрично относительно точки x = 0.Можно показать, что весовые коэффициенты qj (1.2.11) квадратурнойформулы Гаусса положительны (см., например, [3]). Кроме того, в симметричных относительно точки x = 0 корнях полинома Лежандра xj = −xn−(j−1)весовые коэффициенты, соответствующие этим узлам, совпадают при любомn: qj = qn−(j−1) .17Глава 1.
Численные методы вычисления определенного интегралаПриведем значения корней xi и соответствующих им весов qi квадратурных формул Гаусса для n = 1, ..., 5:n=1:n=2:n=3:n=4:n=5:x1 = 0, q1 = 2;p−x1 = x2 = 1/3, q1 = q2 = 1;p−x1 = x3 = 3/5, x2 = 0, q1 = q3 = 5/9, q2 = 8/9;qq√√−x1 = x4 = (15 + 2 30)/35, −x2 = x3 = (15 − 2 30)/35,√√q1 = q4 = (18 − 30)/36,q2 = q3 = (18 + 30)/36;qq√√−x1 = x5 = (35 + 2 70)/63, −x2 = x4 = (35 − 2 70)/63,√√q1 = q5 = (322 − 13 70)/900, q2 = q4 = (322 + 13 70)/900,x3 = 0, q3 = 128/225.(1.2.13)Численные значения узлов xi и весов qi (1.2.13) с десятью десятичнымизнаками после запятой приведены в табл.
1.1. На оценке погрешности квадратурных формул Гаусса останавливаться не будем (см., например, [5]).1.2.7. Правило Рунге практической оценки погрешностиПри выводе формулы средних прямоугольников предполагалось, что f ∈C 2 [a, b]. Погрешность этой формулы, выражающаяся через вторую производную f 00 (x), есть величина O(h2 ). Если подынтегральная функция имеет производные более старших порядков, то можно получить более содержательнуюоценку погрешности.Если f ∈ C 4 [a, b], то можно получить следующее выражение для18Глава 1. Численные методы вычисления определенного интегралаI=Rbf (x)dx:aI = Ihnp + c · h2 + O(h4 ),(1.2.14)где Ihnp — значение интеграла, вычисленное по составной формуле среднихпрямоугольников с шагом h (h = (b − a)/n); c — постоянная, не зависящаяRb 001от h, c = 24 f (x)dx.aВеличина ch2 в выражении (1.2.14) называется главной частью погрешности формулы средних прямоугольников.