36. Интерполяционная формула Ньютона.

Рассмотрим случай табличного задания функции одной переменной. С математической точки зрения это означает, что значения функции f(x) на промежутке [a, b] заданы в узлах одномерной сетки

hx = {xi / xi = xi – 1 + hi, hi > 0, i = 1, 2, 3, …, n; x0 = a, xn = b},         (3.1.1)

где xi, i = 0, 1, 2, …, n – узлы одномерной сетки hx; hi, i = 1, 2, 3, …, n – шаг сетки hx. Если hi = h " i = 1, 2, 3, …, n, сетка hx равномерная.

 

В случае табличного задания функции на отрезке [a, b] часто возникает проблема восстановления значений функции f(x) в произвольной точке x Î (a, b) между узлами сетки hx. Такую задачу называют задачей интерполирования (или интерполяции) значений функции в точку x Î (a, b).

 

В более общей постановке ставится задача замены таблично заданной на отрезке [a, b] функции f(x) некоторой функцией j (x) из заданного класса функций (например, полиномов той или иной степени, экспоненциальных функций и т. д.). При этом выбор функции j (x) из заданного класса осуществляется таким образом, чтобы ее значения были в наибольшей степени согласованы в смысле некоторого критерия с табличными данными. Такая задача называется задачей аппроксимации или приближения таблично заданных функций. Формально задачи интерполяции таблично заданной функции решаются посредством аппроксимации ее функциями заданного класса.

 

Приведем формальную постановку задачи полиномиальной интерполяции таблично заданной на отрезке [a, b] функции f(x). Пусть на отрезке [a, b] задана одномерная сетка

 

hx = {xi / xi = x– 1 + hi, hi > 0, i = 1, 2, 3, …, n; x0 = a, xn = b},

в узлах xi которой заданы yi = f(xi), i = 0, 1, 2, …, n – соответствующие значения функции f(x).

 

Функция j (x) называется интерполирующей или интерполяционной для f(x) на [a, b], если ее значения j (x0), j (x1), …, j (xn) в узлах сетки
x0, x1, x2,…, xn совпадают с заданными значениями yi функции f(xi), i = 0, 1, 2, …, n. Геометрически это означает, что график функции j (x) проходит так, что по крайней мере в точке (n + 1) (в точках, соответствующих узлам сетки hx) на отрезке [a, b] он пересекает или касается истинного графика функции f(x).

 

Интерполяционный полином в форме Ньютона для x Î (a, b) может быть записан в виде

 

Pn(x) = f(x0) + f(x0, x1)×(xx0) + f(x0, x1, x2)×(xx0)(xx1) + …

…+ f(x0, x1, x2, …, xn – 1, xn)×(xx0)(xx1) … (xxn - 2)(xxn - 1)    (3.3.1)

или в более компактной форме

 

,                                 (3.3.2)

где f(x0, x1,…, xk)разделенная разность k-го порядка.

 

 

* Разделенная разность — обобщение понятия производной для дискретного набора точек. Разделенная разность нулевого порядка функции f(x) — сама функция f(x). Разделённая разность порядка n определяется через разделённую разность порядка n − 1 по формуле

f(x_0;\;x_1;\;\ldots;\;x_n)=\frac{f(x_1;\;\ldots;\;x_n)-f(x_0;\;\ldots;\;x_{n-1})}{x_n-x_0}.

Для разделённой разности также верна формула

f(x_0;\;x_1;\;\ldots;\;x_n)=\sum_{j=0}^n\frac{f(x_j)}{\prod\limits_{i=0\atop i\neq j}^n(x_j-x_i)}.

Из этой формулы следует, что разделённая разность является симметрической функцией своих аргументов (то есть при любой их перестановке не меняется), 

 

 

Разделенные разности для табличных данных рассчитываются следующим образом:

·          разделенные разности 1-го порядка по формуле

·          разделенные разности 2-го порядка через разделенные разности 1-го порядка:

 

В общем случае разделенные разности k-го порядка определяются
через разделенные разности (k – 1)-го порядка:

 

Учитывая связь разделенных разностей с производными таблично заданной функции f(x), интерполяционный полином в форме Ньютона можно считать разностной аппроксимацией формулы Тейлора для функции f(x).

 

Нетрудно показать, что для интерполяционного полинома Pn(x)
в форме Ньютона выполняются условия интерполяции по узлам сетки:

Pn(xk) = f(xk), k = 0, 2, …, n.

 

Кроме того, можно убедиться в том, что интерполяционный полином в форме Ньютона только по форме отличается от интерполяционного полинома в форме Лагранжа, представляя собой на заданной сетке один и тот же интерполяционный полином. Поскольку, как было отмечено ранее, интерполяционный полином единственен, то последнее утверждение является вполне естественным и должно выполняться для любых форм представления этого полинома, если сетка и значения функции не изменяются.

 

В соответствии с соотношением (3.3.1) процесс интерполяции таблично заданной на отрезке [a, b] функции в заданную точку x* Î (a, b)
с помощью интерполяционного полинома в форме Ньютона можно представить в виде следующего алгоритма:

 

0. На отрезке [a, b] задать одномерную сетку

hx = {xi / xi = xi – 1 + hi, hi > 0, i = 1, 2, 3, …, n; x0 = a, xn = b}

и значения yi = f(xi) в узлах сетки xi, i = 0, 1, 2, …, n (n ≥ 1).

Задать x* Î (a, b).

1. Положить S: = f(x0); D: = 1; k: = 1.

2. Вычислить разделенные разности k-го порядка

f(xi, xi + 1, xi + 2, …, xi + k) для всех i = 0, 1, 2, …, nk;

D: = D×(x*xk – 1).

3. Вычислить r = f(x0, x1, x2, …, xk)×D.

4. Вычислить S: = S + r.

5. Положить k: = k + 1.

6. Если k £ n, перейти к выполнению п. 2.

7. Процесс завершен. S = Pn(x*) – значение интерполяционного полинома в точке x*.

 

Отметим некоторые характерные особенности использования интерполяционного полинома в форме Ньютона для интерполяции табличных данных.

 

1. Вычисления разделенных разностей в п. 3 приведенного алгоритма удобно производить, записывая их в таблицу, аналогичную табл. 3.3.1 для сетки с девятью узлами.

Таблица 3.3.1

№ узла

xi

f(xi)

f(xi, xi+1)

f(xi, xi+1, xi+2)

f(xi, xi+1, xi+2, xi+3)

f(xi, xi+1, xi+2, …, xi+n)

0

x0

f(x0)

 

 

 

 

 

 

 

f(x0, x1)

 

 

 

1

x1

f(x1)

 

f(x0, x1, x2)

 

 

 

 

 

f(x1, x2)

 

f(x0, x1, x2, x3)

 

2

x2

f(x2)

 

f(x1, x2, x3)

 

 

 

 

 

f(x2, x3)

 

f(x1, x2, x3, x4)

 

3

x3

f(x3)

 

f(x2, x3, x4)

 

 

 

 

 

f(x3, x4)

 

f(x2, x3, x4, x5)

 

4

x4

f(x4)

 

f(x3, x4, x5)

 

f(x0, x1, x2,…, x8)

 

 

 

f(x4, x5)

 

f(x3, x4, x5, x6)

 

5

x5

f(x5)

 

f(x4, x5, x6)

 

 

 

 

 

f(x5, x6)

 

f(x4, x5, x6, x7)

 

6

x6

f(x6)

 

f(x5, x6, x7)

 

 

 

 

 

f(x6, x7)

 

f(x5, x6, x7, x8)

 

7

x7

f(x7)

 

f(x6, x7, x8)

 

 

 

 

 

f(x7, x8)

 

 

 

8

x8

f(x8)

 

 

 

 

Разделенные разности образуют в таблице своеобразный треугольник. При расчете значения Pn(x*) используются разделенные разности, находящиеся как бы на верхней стороне этого треугольника (они выделены в таблице жирным шрифтом).

 

2. Использование интерполяционного полинома в форме Ньютона оказывается удобным, например, когда появляются дополнительные узлы. Так, если добавлен один новый узел, то в формуле (3.3.1) для интерполяционного полинома в форме Ньютона добавится еще одно слагаемое, которое нужно прибавить к уже рассчитанной сумме предыдущих слагаемых (они, очевидно, не изменяются). Для расчета значения дополнительного слагаемого в нижней части таблицы для конечных разностей добавляется еще одна строка, соответствующая этому узлу. Ранее рассчитанные разности пересчитывать не нужно, необходимо лишь определить дополнительно значения тех разделенных разностей, которые будут находиться на нижней стороне нового "треугольника" разделенных разностей.

 

3. Формулу (3.3.1) для интерполяционного полинома в форме Ньютона можно записать для любого способа упорядочивания узлов. Так, например, при упорядочивании узлов в обратном порядке (xn, xn – 1, xn – 2, …, x2, x0) формула (3.3.1) для представления интерполяционного полинома в форме Ньютона примет вид

 

Pn(x) = f(xn) + f(xn, xn – 1)×(xxn) + f(xn, xn – 1, xn – 2)×(xxn)(xxn – 1)+…

…+ f(xn, xn – 1, xn – 2,… x1, x0)×(xxn)(xxn – 1)… (xx2)(xx1).     (3.3.3)

 

Если интерполяционный полином в форме Ньютона Pn(x) определяется по формуле (3.3.1), его называют интерполяционным полиномом для интерполирования вперед (или для интерполирования в начале таблицы). Этот полином целесообразно использовать, когда точка x* находится в начале таблицы.

 

Если Pn(x) определяется по формуле (3.3.3), его называют интерполяционным полиномом для интерполирования назад (или для интерполирования в конце таблицы). Его целесообразно использовать, когда точка x* находится в нижней части таблицы.

 

 

Существенным недостатком полиномиальной интерполяции является то, что при большом числе узлов сетки степень полинома Pn(x), определяемого на отрезке [a, b], будет большой. Известно, что полиномы с большими степенями крайне неустойчиво ведут себя в промежутках между узлами сетки hx, часто принимая там неоправданно большие по абсолютной величине значения, заведомо не соответствующие возможным значениям функции f(x). Поэтому на практике не рекомендуется строить интерполяционные полиномы выше 5-6-й степени.