Интерполяция функций: полиномы Лагранжа, Ньютона и сплайны

Термометр писал температуру каждый час: в 12:00 — 21.4°, в 13:00 — 22.1°, в 14:00 — 22.6°. А сколько было в 13:30? Таблица молчит между точками. Единственный честный источник — сама таблица, и интерполяция предлагает контракт: проведём через все точки одну гладкую кривую и спросим значение у неё. Вопрос лишь в том, какую кривую выбрать — и почему наивный выбор иногда разваливается на осцилляции.

Постановка: по таблице восстановить функцию#

Даны $n+1$ точек $(x_0,y_0),\dots ,(x_n,y_n)$ с различными узлами $x_i$. Требуется функция $P$, проходящая через каждую: $P(x_i)=y_i$ — и по ней оцениваются промежуточные значения. Классический ответ — интерполяционный полином степени не выше $n$: у него ровно $n+1$ коэффициентов, по одному условию на каждую точку, система квадратная.

Где это живёт? Калибровочная кривая датчика: пять замеров по эталонам — и прибор умеет переводить напряжение в температуру между точками. Таблицы логарифмов, которыми пользовались до калькуляторов, снабжались колонками поправок — это чистая интерполяция. Анимация и компьютерная графика: ключевые кадры заданы, между ними кривую достраивает интерполяция. Везде одна и та же математика — отличается только цена ошибки: где-то испорченный график, а где-то траектория манипулятора.

Почему через n+1 точку проходит ровно один полином#

Единственность доказывается в три строки. Пусть нашлись два полинома $P$ и $Q$ степени не выше $n$, проходящие через все точки. Их разность $R=P-Q$ — тоже полином степени не выше $n$, и она обнуляется во всех $n+1$ узлах. Но ненулевой полином степени $n$ имеет не более $n$ корней. Противоречие — значит, $R\equiv 0$, и полином один. Существование докажет формула Лагранжа: она предъявит кандидата явно.

Полином Лагранжа: собираем базис#

Лагранж не решает систему, а собирает ответ из деталей. Для каждого узла строится базисный полином $l_i(x)$ — тот, который равен единице в своём узле и нулю во всех остальных. Тогда сумма $\sum y_i{l}_i(x)$ автоматически проходит через всю таблицу:

Знаменатель ${\prod }_{j\ne i}(x_i-x_j)$ не равен нулю, потому что узлы различны, — так что базис всегда построится. Разберём полный пример с тремя точками: $(0,2)$, $(1,4)$, $(2,3)$.

1. 1
   Базис первого узла: $l_0(x)=\frac{(x-1)(x-2)}{(0-1)(0-2)}=\frac{(x-1)(x-2)}{2}$
2. 2
   Базис второго: $l_1(x)=\frac{x(x-2)}{(1-0)(1-2)}=-x(x-2)$; третьего: $l_2(x)=\frac{x(x-1)}{(2-0)(2-1)}=\frac{x(x-1)}{2}$
3. 3
   Собираем: $L(x)=2\cdot \frac{(x-1)(x-2)}{2}-4x(x-2)+3\cdot \frac{x(x-1)}{2}$
4. 4
   Раскрываем скобки: $L(x)=-1.5x^2+3.5x+2$
5. 5
   Проверка в узлах: $L(0)=2$, $L(1)=4$, $L(2)=3$ — все три условия выполнены

Проверка в узлах — обязательный ритуал: если хоть одно значение не сошлось, где-то потерян знак. Полиномом можно пользоваться и дальше: $L(0.5)=-0.375+1.75+2=3.375$ — оценка температуры в полдень между замерами.

Полином Ньютона: таблица разделённых разностей#

Ньютон идёт другим путём: строит полином последовательно, как лестницу, на ступеньках которой стоят разделённые разности. Разность первого порядка: $f[x_i,x_{i+1}]=\frac{y_{i+1}-y_i}{x_{i+1}-x_i}$ — наклон секущей. Разности второго порядка считаются из разностей первого, и так до вершины. Для нашей таблицы:

Подставляем: $N(x)=2+2x-1.5x(x-1)$. Раскройте скобки — получится тот же $-1.5x^2+3.5x+2$. Полиномы Лагранжа и Ньютона — это одно и то же решение в разных одеждах; единственность гарантирует совпадение. Зачем тогда Ньютон? Смотрите, что происходит при добавлении точки $(3,1)$: считаем лишь новые разности — $f[2,3]=-2$, $f[1,2,3]=-0.5$, $f[0,1,2,3]=\frac{1}{3}$ — и дописываем один член $+\frac{1}{3}x(x-1)(x-2)$. Всё, что было посчитано раньше, остаётся в силе. Лагранжу пришлось бы пересобирать весь базис заново.

Проблема Рунге: почему равномерная сетка подставляет#

Интуиция подсказывает: больше точек — точнее полином. Рунге нашёл функцию, на которой интуиция ломается: $f(x)=\frac{1}{1+25x^2}$ на $[-1;1]$. Возьмём 11 узлов с шагом $0.2$ и построим интерполяционный полином десятой степени:

Максимум ошибки достигается у краёв — около $x=\pm 0.94$ — и составляет $1.92$, хотя сама функция всюду меньше единицы. Причина в поведении произведения $\omega (x)={\prod }_{j=0}^n(x-x_j)$, которым оценивается погрешность интерполяции: на равномерной сетке у краёв это произведение разрастается, и чем больше узлов, тем хуже. Добавление узлов не лечит, а усугубляет: степень растёт, осцилляции у краёв становятся выше. Лечатся два пути: распределить узлы кучнее у краёв — узлы Чебышёва $x_k=\cos \frac{(2k+1)\pi }{2n+2}$ гасят рост $\omega$ — или отказаться от одного глобального полинома.

Как это выглядит на нашем учебном примере с тремя узлами $0,1,2$? Ошибка в точке $x$ пропорциональна $\frac{|f^{\prime \prime \prime }(\xi )|}{6}|x(x-1)(x-2)|$. В точке $x=0.5$ произведение $\omega (0.5)=0.5\cdot (-0.5)\cdot (-1.5)=0.375$ — крошечное, потому и полином второй степени там честен. Вся геометрия узлов живёт в $\omega$: на равномерной сетке у краёв оно в разы больше, чем в центре, а с ростом числа узлов — растёт лавинообразно. Маленькое $\omega$ — щит интерполяции; узлы Чебышёва минимизируют его максимум по отрезку — в этом и есть их секрет.

Кубические сплайны: кусочные кубики с гладким стыком#

Сплайн — гибкая рейка чертёжника: между узлами функция кусочно-кубическая, но стыки согласованы так, что первая и вторая производные непрерывны — гладкость $C^2$. Считаем условия: на $n$ отрезках по 4 коэффициента кубики — всего $4n$ неизвестных. Проход через $n+1$ узлов даёт $n+1$ условий, непрерывность $S^{\prime }$ и $S^{\prime \prime }$ во внутренних стыках — по $n-1$, плюс два условия на краях (например, естественный сплайн $S^{\prime \prime }=0$ на концах). Итого $n+1+2(n-1)+2=4n$ — система решается, причём трёхдиагональным алгоритмом за линейное время.

Почему сплайн не осциллирует? Кубик на каждом отрезке видит только своих четырёх соседей — ошибка не успевает накопиться и расползтись по всему отрезку, как у глобального полинома. Числа для той же функции Рунге: естественный кубический сплайн по тем же 11 узлам даёт максимум ошибки $0.022$ против $1.92$ у полинома десятой степени — разница почти в девяносто раз. В точке $x=0.95$ сплайн возвращает $0.0429$ против истинных $0.0424$: три верных значащих цифры там, где полином выдавал $1.92$.

Слово «сплайн» пришло из чертёжной практики: гибкая рейка, прижатая гвоздиками к опорным точкам, изгибается так, как требует минимум энергии — кусочные кубики с непрерывной кривизной. Отсюда и слава сплайнов в инженерии: кузова автомобилей, обводы кораблей, траектории станков — везде, где нужна гладкость второй производной, сплайн — рабочий инструмент, а не только учебная задача.

Когда интерполяция уместна, а когда нужен МНК#

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

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