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

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

Форма Ньютона записывает интерполяционный полином как $P(x) = c_0 + c_1(x-x_0) + c_2(x-x_0)(x-x_1) + \ldots$, где коэффициенты $c_k$ — разделённые разности.

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

Полином через узлы — тот же самый, что у Лагранжа (он единственный), но форма записи другая. Ньютон складывает его из слагаемых, каждое из которых добавляет один узел и не портит совпадение в предыдущих. Первое слагаемое — константа $c_0 = y_0$. Второе — $c_1(x-x_0)$: оно зануляется в $x_0$ (там уже всё верно) и подправляет значение в $x_1$. Третье — $c_2(x-x_0)(x-x_1)$: зануляется в $x_0$ и $x_1$, подправляет $x_2$. И так далее. Главный плюс: добавили узел — дописали одно слагаемое, ничего не пересчитывая.

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

Коэффициенты $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_2)(x-x_1) + c_1)(x-x_0) + c_0$. Разделённые разности симметричны относительно перестановки узлов, поэтому порядок узлов на итоговый полином не влияет (хотя на численную устойчивость промежуточных вычислений — немного влияет). Связь с производными прямая: $k$-я разделённая разность равна $f^{(k)}(\xi)/k!$ в некоторой точке $\xi$ — поэтому форма Ньютона так похожа на ряд Тейлора, только «по узлам».

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

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

Итоги

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