1626435587-51311eae4652e8ad616b5bdef025cbb3 (844239), страница 15
Текст из файла (страница 15)
Действуяаналогично дальше, можно показать, что разделённая разность порядка является полиномом степени − при = 1, 2, . . . , . Для разделённой разности порядка + 1 имеем: (, 0 , 1 , . . . , ) = (, 0 , . . . , −1 ) − (0 , . . . , ). − (52)При = числитель дроби зануляется в силу симметрии разделённыхразностей по перестановкам аргументов. С другой стороны, выражение74в числителе дроби является полиномом степени 0, т. е. константой.Следовательно, числитель дроби тождественно равен нулю ∀.Последовательно подставляя друг в друга выражения (50), (51), .
. . ,(52) и учитывая на последнем шаге, что (, 0 , 1 , . . . , −1 ) − (0 , 1 , . . . , −1 , ) = 0,приходим к полиному, записанному по () = (0 )...схеме Горнера :+ ( − 0 ) · [ (0 , 1 ) ++ ( − 1 ) · [ (0 , 1 , 2 ) + . . .+ ( − −2 ) · [ (0 , . . . , −1 ) ++ ( − −1 ) · (0 , 1 , .
. . , )] . . . ]] .(53)Или, раскрывая все квадратные скобки в (53), можно записать формулу Ньютона в виде суммы произведений: () = (0 ) +∑︁ (0 , 1 , . . . , )−1∏︁( − ).(54)=0=1Для вычисления полученных выражений необходимо заметить, чторазделённые разности для полиномов (0 , 1 , . . . , ) совпадают с разделёнными разностями для интерполируемой функции (), поскольку ( ) = ( ) ∀ = 0, . . . , . Следовательно, формулы (53) и (54)после замены (0 , 1 , . . .
, ) на (0 , 1 , . . . , ) дают решение задачи о построении интерполяционного полинома в форме Ньютона. Прификсированном числе узлов сетки для вычислений более удобна схема Горнера (53). Напротив, формула Ньютона (54) предпочтительнее вслучае, когда необходимо контролировать точность и добавлять новыеузлы в расчёт. При этом полезно помнить о том, что порядок, в котором перенумерованы узлы интерполяции 0 , 1 , .
. . , не играет роли всилу симметрии разделённых разностей по перестановкам аргументов.Для вычисления разделённых разностей, входящих в выражения(53) и (54), удобно записывать коэффициенты в виде таблицы (в виде соответствующего двумерного массива при вычислении на ЭВМ).Для последующего вычисления значений интерполяционных полиномов необходимо сохранять в памяти ЭВМ только верхнюю косую строку табл. 2.Число операций, необходимых для вычисления значения интерполяционного полинома, записаного по схеме Горнера (53), равно 3:на каждом из шагов необходимо выполнить два сложения и одно75умножение. Обратим внимание, что число операций в формуле Ньютона растёт линейно по , а не квадратично, как в формуле Лагранжа(48).
Из сравнения полученных ответов для числа операций видно, чтоформула Ньютона более экономична по сравнению с формулой Лагранжа для любой степени интерполяционного полинома.Таблица 20123401234(0 , 1 )(1 , 2 )(2 , 3 )(3 , 4 )(0 , 1 , 2 )(1 , 2 , 3 )(2 , 3 , 4 )(0 , 1 , 2 , 3 )(1 , 2 , 3 , 4 )(0 , 1 , . . . , 4 )Для построения табл. 2 (нахождения коэффициентов полинома)необходимо вычислить + ( − 1) + . . . + 1 = 12 ( + 1) разделённыхразностей; для вычисления каждой из них требуется два вычитания иодно деление, итого 23 ( + 1) арифметических операций. Как и в случае формулы Лагранжа, число «подготовительных» операций растёткак (2 ), хотя численный коэффициент для формулы Ньютона получился приблизительно на четверть меньше.
Заметим, однако, что этиоперации выполняются однократно, так что временем на их выполнение в большинстве задач можно пренебречь.Заметим также, что нахождение полиномиальных коэффициентовв формуле Ньютона (53) возможно и без введения понятия разделённых разностей. Действительно, запишем интерполяционный полином Ньютона по схеме Горнера с неопределёнными коэффициентами0 , 1 , . . . , : () = 0 + ( − 0 ) [1 + . . . + ( − −2 ) [−1 + ( − −1 ) · ] .
. .] .Требуя совпадения значений полинмома с интерполируемой функциейв узлах сетки ( ) = ( ) ≡ , получим систему линейных алгебраических уравнений на коэффициенты 0 , 1 , . . . , :012===...=00 + (1 − 0 )10 + (2 − 0 ) [1 + (2 − 1 )2 ]0 + ( − 0 ) [1 + . . . + ( − −2 ) [−1 ++ ( − −1 ) ] . . .]Из первого уравнения имеем 0 = 0 ; подставляя 0 во второе, находим01 = 11 −−0 ; далее подставляем 0 , 1 в третье уравнение и т. д.765.3. Погрешность интерполяцииОценим погрешность интерполяции полиномами в случае, когда интерполируемая функция имеет + 1 непрерывную производную. Поскольку погрешность интерполяции обращается в ноль во всех узлахсетки 0 , 1 , .
. . , , её можно представить в следующем виде:() − () = () · (),() =∏︁( − ).(55)=0Введём вспомогательную функцию() ≡ () − () − ()(),(56)где играет роль параметра. Очевидно, что функция () обращаетсяв ноль в каждом узле сетки 0 , 1 , . . . , и в точке = . Посколькупо нашему предположению () имеет n+1 непрерывную производную,() также должна быть +1 раз непрерывно дифференцируемой, причём для ( + 1)-й производной в силу определения () можно записать (+1) () = (+1) () − ( + 1)! ().(57)Из курса анализа известно, что между двумя нулями гладкой функции лежит ноль её производной. Последовательно применяя это правило, можно сделать вывод, что в промежутке между крайними из + 2нулей функции () лежит ноль её ( + 1)-й производной.
Обозначимточку, в которой зануляется производная (+1) () за * = * () (онабудет являться некоторой функцией от параметра ). Приравнивая (57)к нулю, получим (+1) (* ())() =.(58)( + 1)!Поскольку параметр изначально был выбран произвольно, равенство(58) позволяет сделать оценку сверху для погрешности полиномиальной интерполяции (55):|() − ()| ≤max | (+1) ()|· |()|.( + 1)!а(59)Для примера на рис. 6 ( ) показан график зависимости разности() − () при интерполяции функции () = sin() полиномами степени = 5 на промежутке от 0 до /2, а также зависимость максимальной ошибки от степени (рис. 6 ( )).77б(б) 100ω(x)r(x)ω(x)/6!210-3НевязкаНевязка x 105(а)1010-610-910-12-110-1500.5101.5481216NXРис.
6. Зависимость погрешности от(а) и степени полинома(б)На графике видно, что погрешность зануляется в шести точках0 , 1 , . . . , 5 . Ближе к краям отрезка [0, /2], на котором производится интерполяция, ошибка больше, чем в середине. За пределами отрезка ошибка резко возрастает — погрешность экстраполяции намногопревышает погрешность интерполяции. В целом погрешность близка кмажорантной оценке (59) (показана на рис. 6 ( ) серой линией).На рис.
6 ( ) приведена также зависимость максимальной (по прификсированном ) ошибки как функции степени полинома (по оси использован логарифмический масштаб). Видно, что ошибка быстроубывает до = 13, после чего вновь начинает возрастать из-за погрешности вычисления разделённых разностей (см. п. 5.4).аб5.4. Численное дифференцированиеПокажем, как разделённые разности, введённые на с.
73, могутбыть использованы для аппроксимации значений производных функции. Пусть () — достаточное число раз гладкая функция, для которойсправедливо разложение() =∑︁=0 () ()(︀)︀( − )+ ( − )+1!(60)Вычислим разделённую разность -го порядка (0 , 1 , . . . , ), переписав (60) в виде () = () + (ℎ+1 ), где () — полином степени, ℎ ≤ max | − |/ — характерный шаг сетки {0 , .
. . , }:(0 , . . . , ) = (0 , . . . , ) + (ℎ).78(61)Для вычисления разделённой разности от величины (ℎ+1 ) мы воспользовались формулой (49) на с. 74, получив (ℎ).Входящая в (61) разделённая разность -го порядка для полиномастепени равна -й производной ( / ) этого полинома, делённойна ! (в этом легко убедиться, продифференцировав раз выражение(54) на с. 75). Заметим, что тот факт, что выражение (54) было получено нами для интерполяционного полинома, абсолютно несущественен:в силу доказанной на с.
72 теоремы о существовании и единственности,любой полином можно считать интерполяционным по отношению ксамому себе. Учитывая, что ( / ) = ( / ) в силу определения (60), получаем из (61) окончательно: () () = ! (0 , . . . , ) + (ℎ).(62)Выбор точки достаточно произволен — она может совпадать с любымиз узлов сетки или находиться на расстоянии (ℎ) от него. Порядок точности аппроксимации можно повысить, используя формулы сбо́льшим числом точек.
(В общем случае произвольной неравномернойсетки порядок точности по ℎ равен разности числа узлов и порядкапроизводной [2, гл. 3].)Разделённую разность порядка , вычисленную на равномерной сетке и умноженную на !, называют. Использованиеравномерной сетки обычно позволяет упростить формулы и повыситьпорядок точности аппроксимации. Для примера рассмотрим конечнуюразность второго порядка:конечной разностью2(ℎ, 0, −ℎ) =(−ℎ) − 2(0) + (ℎ)ℎ2= ′′ (0) + (iv) (0) + . .
.2ℎ12(63)Конечная разность второго порядка (63) с точностью до (ℎ2 ) аппроксимирует значение второй производной функции в среднем узле. Аналогично, конечная разность третьего порядка)︁1 (︁3! (0, ℎ, 2ℎ, 3ℎ) = 3 (3ℎ) − 3(2ℎ) + 3(ℎ) − (0)ℎаппроксимирует ′′′ (3ℎ/2) с точностью ≈ ℎ2 (v) /8 = (ℎ2 ) и т. д.Заметим, что конечные разности нечётного порядка аппроксимируют значения соответствующих производных со вторым порядком точности в средней точке между центральными узлами сетки, тогда какаппроксимация в узлах имеет лишь первый порядок по ℎ. Как следствие, для аппроксимации производных нечётного порядка в узлах ,конечные разности выгодно вычислять на сетке с удвоенным шагом 2ℎ.79Хотя это и увеличивает численный коэффициент в остаточном члене в2 раз, но обеспечивает более высокий порядок точности по ℎ.Несмотря на то, что из выражения (62) следует (0 , .
. . , ) → ()в пределе max | − | → 0, вычисление производных, особенно высокихпорядков, может представлять существенные сложности. Действительно, пусть требуется аппроксимировать величину () , вычисляя конечную разность ! (0 , . . . , ), = 0 + ℎ. Изменим входные данные —значения ≡ ( ) — на малую величину · exp( ), || ≪ | |. Приэтом конечная разность, вычисленная по изменённым данным, очевидно, будет аппроксимировать величину ( / )[() + exp()] = () + () exp().