Численные методы. Ионкин (2012) (неоффициальные) (косяки есть) (1160437), страница 8
Текст из файла (страница 8)
. . , xk ) =0 (x )wii=0где w(x) = (x − x0 )(x − x1 ) . . . (x − xn ). Будем обозначать через wα,β (x) :wα,β (x) = (x − xα )(x − xα+1 ) . . . (x − xβ )0Тогда, используя равенство w0 (xi ) = w0,k(xi ) =kQ(xi − xj ) можно записатьj=0,j6=iусловие утверждения следующем образом:f (x0 , x1 , . . . , xk ) =kXf (xi )0w0,k(xi )i=0(1)Доказательство: Доказательство проведем методом полной математической индукf (x1 ) − f (x0 )f (x0 )f (x1 )ции.
По определению f (x0 , x1 ) =. Тогда f (x0 , x1 ) =+,x1 − x0x0 − x1 x1 − x0то есть для k = 1 формула справедлива. Предположим, что она справедлива и дляk = l:l+1lXXf (xi )f (xi ), f (x1 , . . . , xl , xl+1 ) =(2)f (x0 , x1 , . . . , xl ) =00w0,l (xi )w1,l+1 (xi )i=1i=0По определению составим l + 1 разделенную разность и представим её в виде трехслагаемых:f (x0 )f (xl+1 )++00(x0 − xl+1 )w0,l (x0 ) (xl+1 − x0 )w1,l+1(xl+1 )!lXf (xi )11+− 00x − x0 w1,l+1(xi ) w0,l(xi )i=1 l+1f (x0 , x1 , .
. . , xl+1 ) =0000Очевидно, что: (xl+1 − x0 )w1,l+1(xl+1 ) = w0,l+1(xl+1 ) и (x0 − xl+1 )w0,l(x0 ) = w0,l+1(x0 ).Далее рассмотрим выражение , стоящее под знаком суммы:!111− 00xl+1 − x0 w1,l+1(xi ) w0,l(xi )54Домножая числитель и знаменатель первой дроби, стоящий в скобках на xi − x0 , авторой — на xi − xl+1 , получим:!1xi − x0xi − xl+11−=000xl+1 − x0 w1,l+1(xi )(xi − x0 ) w0,l(xi )(xi − xl+1 )w0,l+1(xi )0000Очевидно, что w1,l+1(xl+1 )(xi − x0 ) = w0,l+1(xi ), w0,l(xi )(xi − xl+1 ) = w0,l+1(xi ).Подставляя полученные выражения в исходную сумму, получаем формулу (1) дляk = l + 1:l+1Xf (xi )f (x0 , x1 , . . .
, xl+1 ) =0w0,l+1 (xi )i=0Выразим значение функции в k−м узле через f0 и разделенные разности доk−го порядка:при k = 1 : f (x0 , x1 ) =f (x1 ) − f (x0 )f (x0 )f (x1 )=+x1 − x0x 0 − x1 x1 − x0Домножив обе части на x1 − x0 получим:(x1 − x0 )f (x0 , x1 ) = f (x1 ) − f (x0 ) и f (x1 ) = f (x0 ) + (x1 − x0 )f (x0 , x1 )при k = 2 : f (x0 , x1 , x2 ) =(3)f (x1 )f (x2 )f (x0 )++(x0 − x1 )(x0 − x2 ) (x1 − x0 )(x1 − x2 ) (x2 − x1 )(x2 − x0 )Следовательно:(x2 − x1 )(x2 − x0 )f (x0 , x1 , x2 ) =f (x0 )(x2 − x1 ) f (x1 )(x2 − x0 )++ f (x2 )x 1 − x0x0 − x1Подставим ранее полученное для f (x1 ) и объединим полученные слагаемые:f (x2 ) = f (x0 ) + f (x0 , x1 )(x2 − x0 ) + (x2 − x0 )(x2 − x1 )f (x0 , x1 , x2 )(4)Тогда можно сказать, чтоf (xk ) = f (x0 ) + (xk − x0 )f (x0 , x1 ) + .
. . + (xk − x0 ) . . . (xk − xk−1 )f (x0 , x1 , . . . , xk ) (5)Замечание. Безусловно, формулу (5) можно аккуратно доказать методом полнойматематической индукции. Мы показали лишь переход от k = 1 к k = 2.§4 Интерполяционная формула НьютонаПолучим в явном виде интерполяционный полином Ньютона Nn (x) для функции f (x) по узлам {xi }n0 . Для этого воспользуемся формулой (5). Полином Nn (x)получается из (5) заменой xn на x:Nn (x) = f (x0 ) + (x − x0 )f (x0 , x1 ) + . .
. + (x − x0 ) . . . (x − xn−1 )f (x0 , x1 , . . . , xn )(6)55По определению интерполяционного полинома нужно показать,что Nn (xi ) =f (xi ) ∀i = 0, n. Ясно, чтоNn (xi ) = f (x0 ) + (xi − x0 )f (x0 , x1 ) + . . . + (xi − x0 ) . . . (xi − xi−1 )f (x0 , x1 , . . . , xi )По формуле (5) это равно f (xi ). А значит полином (6) является интерполяционным полиномом и носит название интерполяционного полинома Ньютона.
Длявычисления погрешности интерполяционного полинома Ньютона можно воспользоваться формулой:f (n+1) (ξ)w(x)ψNn (x) =(n + 1)!Или в более привычной форме, положив Mn+1 = sup |f n+1 (x)| по всем x :|ψNn (x) | 6Mn+1|w(x)|(n + 1)!Замечание. Если узлы фиксированы, то удобен полином Лагранжа, а если фиксирована функция, а количество узлов увеличивается на каждой итерации, то удобенполином Ньютона.§5 Интерполирование с кратными узлами. ПолиномЭрмитаПостановка задачи: Пусть заданы m + 1 узел {xi }m0 и значения функциив этих узлах f (x0 ) .
. . f (xm ). Кроме того, в каждом узле заданы производныеf (a0 −1) (x0 ) . . . f (am −1) (xm ) где ai − кратность для узла xi . Задача заключается в построении полинома n−й степени, значения в узлах которого совпадают со значениямизаданной функции, а его производные - со значением соответствующих производныхзаданной функции:Hn(i) (xk ) = f (i) (xk )(1)Тогда ясно, что a0 +a1 +· · ·+am = n+1.
Если это выполнено, то полином Эрмитаищется и представляется в виде:Hn (x) =m aXk −1Xck,i (x)f (i) (xk )(2)k=0 i=0где ck,i (x)− полином n−й степени от x.Построим полином Эрмита с кратными узлами H3 (x), где узел x1 −кратный (в немзаданы значения f (x1 ), f 0 (x1 )), а узлы x0 , x2 − простые:H3 (x0 ) = f (x0 )H (x ) = f (x )3 11H3 (x2 ) = f (x2 ) 0H3 (x1 ) = f 0 (x1 )56Будем искать его в виде:H3 (x) = c0 (x)f (x0 ) + c1 (x)f (x1 ) + c2 (x)f (x2 ) + b1 (x)f 0 (x1 )Ясно,что:c0 (x0 ) = 1c (x ) = 00 1c0 (x2 ) = 0 0c0 (x1 ) = 0c1 (x0 ) = 0c (x ) = 11 1c1 (x2 ) = 0 0c1 (x1 ) = 0c2 (x0 ) = 0c (x ) = 02 1c2 (x2 ) = 1 0c2 (x1 ) = 0(3)b1 (x0 ) = 0b (x ) = 01 1b1 (x2 ) = 0 0b1 (x1 ) = 1Исходя из этой таблицы, выпишем последовательно все коэффициенты полинома H3 (x).Ищем c0 (x) в виде: c0 (x) = k(x − x2 )(x − x1 )2 , а k найдем из условияc0 (x0 ) = 1. Значит 1 = k(x0 − x2 )(x0 − x1 )2 и :c0 (x) =(x − x2 )(x − x1 )2(x0 − x2 )(x0 − x1 )2Аналогичными рассуждениями можно получить коэффициент c2 (x) :c2 (x) =(x − x0 )(x − x1 )2(x2 − x0 )(x2 − x1 )2b1 (x) будем искать в виде b1 (x) = k1 (x − x0 )(x − x1 )(x − x2 ).
Перепишем b1 (x) в виде:b1 (x) = (x − x1 )[k1 (x − x0 )(x − x2 )], b01 (x) = [. . . ] + (x − x1 )[. . . ]0 . И в точке x1производная будет равна : b01 (x1 ) = k1 (x1 − x0 )(x1 − x2 ) = 1.. Откуда получаем, чтоb1 (x) =(x − x0 )(x − x1 )(x − x2 )(x1 − x0 )(x1 − x2 )Коэффициент c1 (x) будем искать в виде x1 : c1 (x) = (x − x0 )(x − x2 )(ax + b). Так как:(c1 (x1 ) = 1 = (x1 − x0 )(x1 − x2 )(ax1 + b)c01 (x1 ) = 0 = a(x1 − x0 )(x1 − x2 ) + (ax1 + b)(2x − x0 − x2 )то12x1 − x0 − x2b=a=−22(x1 − x0 ) (x1 − x2 )(x1 − x0 )(x1 − x2 )x1 (2x1 − x0 − x2 )1+(x1 − x0 )(x1 − x2 )Окончательный вид c1 (x) :(x − x0 )(x − x2 )c1 (x) =(x1 − x0 )(x1 − x2 )(2x1 − x0 − x2 )(x − x1 )1−(x1 − x0 )(x1 − x2 )Таким образом, показано, что коэффициенты H3 (x) находятся в явном виде иоднозначно.57Оценка погрешности полинома Эрмита H3 (x), построенного поузлам x0 , x1 , x2Введем функциюq(s) = f (s) − H3 (s) − Kw(s),гдеs ∈ [x0 , x2 ],(4)x ∈ [x0 , x2 ], иw(s) = (x − x0 )(x − x1 )2 (x − x2 ) — полином 4-й степени.Постоянную K(x) выбираем из условия:q(x) = 0 = f (x) − H3 (x) − Kw(x)Тогда ясно, что константу K нужно брать равнойK=f (x) − H3 (x)w(x)Функция q(x) имеет 4 нуля на отрезке [x0 , x2 ].
Предположим, что функция f (x)гладкая. Значит к ней можно применить теорему Ролля. У полученной функции q 0 (x)— не менее 3х нулей. Так как узел x1 – кратный, то на данном этапе добавляетсяq 0 (x1 ) = 0, откуда следует, что q 0 (x) имеет не менее 4х нулей, а значит, q 00 (x) – неменее 3х нулей, q 000 (x) – не менее 2х нулей, и тогда q IV (x) – не менее одного нуля.Таким образом существует точка ξ, в которой q IV (ξ) = 0.Тогда продифференцируем 4 раза функцию q(s):q IV (s) = f IV (s) − k·4!Согласно теореме Ролля ∃ξ ∈ [a, b], в которой q IV (ξ) = 0. Если обозначитьM4 =sup |f IV (x)|,x0 6x6x2то получим|f (x) − H3 (x)| = |ψH3 (x)| 6M4M4|w(x)| =|(x − x0 )(x − x1 )2 (x − x2 )|4!4!Получили оценку для ПЭ H3 (x).Используя оценку ΨH3 получим, что ∀n ∈ N :ΨHn (x) =f (n+1) (ξ)(x − x0 )a0 (x − x1 )a1 .
. . (x − xm )am(n + 1)!И если обозначитьMn+1 = sup |f (n+1) (x)|,xто|ΨHn (x)| 6Mn+1|w(x)|(n + 1)!(5)58Полином H3 (x) может быть получен из многочлена L3 (x) с помощью предельного перехода. Пусть есть узлы x0 , x1 , x2 , x3 , где x3 — фиктивный узел. По этим 4музлам можно построить L3 (x)Первое слагаемое(x − x1 )(x − x2 )(x − x3 )f (x0 ) + . . .(x0 − x1 )(x0 − x2 )(x0 − x3 )L3 (x) =При сведении к H3 (x) коэффициенты (x − x3 ) → (x − x1 ) и (x0 − x3 ) → (x0 − x1 ).Задача.
Доказать, чтоlim L3 (x) = H3 (x)x3 →x1Решение: Рассмотри полином Лагранжа для функции f :(x − x1 )(x − x2 )(x − x3 )f (x0 ) + . . .(x0 − x1 )(x0 − x2 )(x0 − x3 )L3 (x) =Тогда(x − x1 )2 (x − x2 )f (x0 ) + · · · = H3 (x)x3 →x1(x0 − x1 )2 (x0 − x2 )Подробные выкладки предлагаем провести читателю.lim L3 (x) =§6 Использование полинома H3(x) для оценки погрешности квадратурной формулы СимпсонаРассмотрим интегралRbf (x)dx. Будем вычислять его с помощью квадратурнойaформулы Симпсона.Разобьем отрезок на n частичных сегментовa 6 x0 < x1 < x2 < · · · < xN 6 bтак, чтоxi − xi−1 = h = const.Тогда, квадратурная формула Симпсона на частичном сегменте имеет вид:Zxif (x)dx =h(fi−1 + 4fi− 1 + fi ),26xi−1где f (xi ) = fi , а fi− 1 — значение функции в середине отрезка [xi−1 , xi ]:2fi− 1 = f (xi−1 + 0.5h)2Если подынтегральная функция имеет видf (x) = a0 + a1 x + a2 x2 + a3 x3 ,(1)59то квадратурная формула Симпсона просто точна (для второй степени — точна попостроению).Докажем, что формула точна для x3 :Zxix3 dx =(x2 − x2i−1 )(x2i + x2i−1 )x4i − x4i−1= i=44xi−1(xi − xi−1 )(xi + xi−1 )(x2i + x2i−1 )h== (xi + xi−1 )(x2i + x2i−1 )44Тогда3 !h 3hx+xii−1(x + 4x3i− 1 + x3i ) =(xi−1 + xi )(x2i−1 − xi xi−1 + x2i ) + 4=26 i−162(xi + xi−1 )(x2i + 2xi xi−1 + x2i−1 )h22=(xi−1 + xi )(xi−1 − xi xi−1 + xi ) +=62 22xi−1 − 2xi xi−1 + 2x2i + x2i + 2xi xi−1 + x2i−1h== (xi + xi−1 )62Zxihhx3 dx.= (xi + xi−1 )3(x2i−1 + x2i ) = (xi + xi−1 )(x2i + x2i−1 ) =124xi−1Это и есть аналитическое выражение интеграла.
Если подынтегральная функциястепени не выше, чем 3, то квадратурная формула Симпсона (КФС) точна. Следовательно, для ПЭ H3 (x) она будет точна.Теперь получим погрешность для КФС: Построим H3 (x) по узлам xi−1 , xi− 1 , xi ,2гдеH3 (xi−1 ) = f (xi−1 );H3 (xi− 1 ) = fi− 1 ;22H3 (xi ) = fi ;0H30 (xi− 1 ) = fi−122f IV (ξ)(x − xi−1 )(x − xi− 1 )2 (x − xi )24!Теперь подынтегральную функцию заменим наΨH3 (x) =f (x) = H3 (x) + ΨH3 (x)60ПолучимZxiZxiZxif (x)dx = {в силу линейности} =H3 (x)dx +ΨH3 (x)dx =xi−1xi−1=h(H3 (xi−1 ) + 4H3 (xi− 1 ) + H3 (xi )) +26xi−1ZxiΨH3 (x)dx =xi−1h(fi−1 + 4fi− 1 + fi ) + Ψi (f )26Тогда ясно, что погрешность КФС на частичном сегменте:ZxiZxihΨH3 (x)dx,Ψi (f ) =f (x)dx − (fi−1 + 4fi− 1 + fi ) =26=xi−1xi−1ОбозначимM4 =sup|f IV (x)|.xi−1 6x6xiТогда можно сказать, чтоM4|Ψi (f )| 64!Zxi(x − xi−1 )(x − xi− 1 )2 (xi − x)dx2xi−1Задача.ZxiПоказать, что(x − xi−1 )(x − xi− 1 )2 (xi − x)dx =2xi−1h5120Решение: Проведем замену в подынтегральном выражении: x = xi−1 + th, t ∈ [0, 1].2212Тогда dx = hdt и x − xi−1 = th, xi − x = h(1 − t), x − xi− 1 = h t −.22Таким образом подставляя эти выражения в интеграл, получаем:Zxi2(x − xi−1 )(xi − x) x − xi− 1 dx =2xi−1h5Z12Z1 15 21h5534(t)(1 − t) t −dt = h2t − t − t + t dt =24412000M4 h5Таким образом |ψi (f )| 6.
Следовательно, погрешность КФС на частич4! 120ном отрезке имеет 5й порядок по h. Запишем погрешность на всем отрезке:Zbf (x)dx −Ψh (f ) =aNXi=1Ψi (f ),61 4M(b−a)h4Ψh (f ) 6,2180 hN = b − a,Таким образом квадратурная формула Симпсона на всем отрезке интегрированияимеет 4й порядок точности по h.§7 Наилучшее среднеквадратичное приближение функцийВведем для начала пространство функций, интегрируемых с квадратом:H = L2 – гильбертово пространство — пространство функций таких, чтоZb∀f ∈ L2f 2 (x)dx < ∞aЧтобы сделать пространство нормированным, введем скалярное произведение:Zb∀f, g ∈ L2f (x)g(x)dx(f, g) =aи нормуZbkf kL2 = 21f 2 (x)dxaПереходим к постановке задачи. Рассмотрим линейно независимые функцииϕ0 (x), ϕ1 (x), .