Разделённые разности Ньютона

Урок про форму Ньютона интерполяционного полинома — гибкую альтернативу Лагранжу, удобную для наращивания узлов.

Форма Ньютона записывает интерполяционный полином как P(x) = c₀ + c₁(x−x₀) + c₂(x−x₀)(x−x₁) + ..., где коэффициенты c_k — разделённые разности.

Идея: строить полином по нарастающей

Полином через узлы — тот же самый, что у Лагранжа (он единственный), но форма записи другая. Ньютон складывает его из слагаемых, каждое из которых добавляет один узел и не портит совпадение в предыдущих. Первое слагаемое — константа c₀ = y₀. Второе — c₁(x−x₀): оно зануляется в x₀ (там уже всё верно) и подправляет значение в x₁. Третье — c₂(x−x₀)(x−x₁): зануляется в x₀ и x₁, подправляет x₂. И так далее. Главный плюс: добавили узел — дописали одно слагаемое, ничего не пересчитывая.

Разделённые разности

Коэффициенты c_k — это разделённые разности, рекурсивный аналог производных по таблице. Разность нулевого порядка — само значение f[x_i] = y_i. Первого — f[x_i,x_{i+1}] = (f[x_{i+1}] − f[x_i])/(x_{i+1} − x_i) (наклон). Каждая следующая строится из двух соседних предыдущего порядка. Удобно считать на месте в одном массиве.

def разделённые_разности(xs, ys):
    n = len(xs)
    c = ys[:]                       # начнём с y-значений
    for j in range(1, n):
        # пересчёт справа налево, чтобы не затереть нужное
        for i in range(n - 1, j - 1, -1):
            c[i] = (c[i] - c[i - 1]) / (xs[i] - xs[i - j])
    return c                        # c[k] — коэффициенты полинома Ньютона

def ньютон_значение(xs, c, x):
    # схема Горнера: P(x) = (...((c_n)(x-x_{n-1}) + c_{n-1})...)
    n = len(xs)
    результат = c[-1]
    for i in range(n - 2, -1, -1):
        результат = результат * (x - xs[i]) + c[i]
    return результат

xs = [0.0, 1.0, 2.0]
ys = [1.0, 3.0, 2.0]
c = разделённые_разности(xs, ys)
print("коэффициенты Ньютона:", [round(v, 4) for v in c])
print(f"P(0.5) = {ньютон_значение(xs, c, 0.5):.6f}")
print(f"P(1.5) = {ньютон_значение(xs, c, 1.5):.6f}")
# совпадает ли с Лагранжем из прошлого урока? Должно — полином один!
print("в узлах:", [round(ньютон_значение(xs, c, t), 4) for t in xs])

Вывод:

коэффициенты Ньютона: [1.0, 2.0, -1.5]
P(0.5) = 2.375000
P(1.5) = 2.875000
в узлах: [1.0, 3.0, 2.0]

Те же P(0.5)=2.375 и P(1.5)=2.875, что у Лагранжа — потому что полином через эти узлы единственный. Различается лишь форма записи и удобство.

Таблица разделённых разностей

На бумаге разделённые разности удобно строить «лесенкой»: в первом столбце — значения y, в каждом следующем — разности соседних элементов предыдущего, делённые на разность соответствующих x. Коэффициенты полинома — верхняя диагональ таблицы.

x_i  f[x_i]   f[.,.]    f[.,.,.]
 0    1
              (3-1)/(1-0)=2
 1    3                 (-1-2)/(2-0)=-1.5
              (2-3)/(2-1)=-1
 2    2
коэффициенты P: c0=1, c1=2, c2=-1.5  (верхняя диагональ)

Как работает под капотом

Вычисление полинома ведут по схеме Горнера — она минимизирует число умножений и численно устойчива: P(x) = ((c₂)(x−x₁) + c₁)(x−x₀) + c₀. Разделённые разности симметричны относительно перестановки узлов, поэтому порядок узлов на итоговый полином не влияет (хотя на численную устойчивость промежуточных вычислений — немного влияет). Связь с производными прямая: k-я разделённая разность равна f^{(k)}(ξ)/k! в некоторой точке ξ — поэтому форма Ньютона так похожа на ряд Тейлора, только «по узлам».

Частые ошибки

Затирать массив при пересчёте. Внутренний цикл по i идёт справа налево, иначе обновлённые значения испортят ещё не использованные.
Считать, что Ньютон точнее Лагранжа. Полином тот же; разница только в форме и инкрементности, а не в точности.
Вычислять полином «в лоб» по степеням. Раскрытие скобок теряет точность; используйте схему Горнера по форме Ньютона.

Итоги

Форма Ньютона: полином как сумма слагаемых, каждое добавляет один узел и не портит предыдущие.
Коэффициенты — разделённые разности, считаются рекурсивно одним проходом по таблице.
Инкрементность: новый узел = одно новое слагаемое без пересчёта; вычисление — по схеме Горнера.
Итоговая кривая совпадает с Лагранжевой — полином через узлы единственный.