Задача интерполяции. Полиномиальная интерполяция

Связь дискретизации и интерполяции

Дискретизация – это замена непрерывной функции на таблицу значений в конечном количестве точек или на каком-либо дискретном множестве. Сеточно-заданная функция - это уже дискретизованная функция.

В данном семинаре нам интересна обратная задача - восстановление непрерывной функции по дискретному набору точек.

Интерполяция — вид аппроксимации, в котором кривая точно проходит через имеющиеся точки данных. Заменяет дискретный ряд значений непрерывной функцией внутри рассматриваемого промежутка.

Основные виды задач интерполяции:

• Нахождение промежуточного значения между двумя точками на кривой

• Построение гладкой кривой, соединяющей точки с известными значениями

• Получение непрерывной функции на основе входных данных, в которой устранено влияние шума

• Нахождение многочлена, хорошим образом приближающим функцию на отрезке (а не только в окрестности точки)

• Экстраполяция дискретно-заданной функции за пределы отрезка, на котором она определена

Интерполяционный многочлен (Лагранжева интерполяция)

Пусть на отрезке \([a, b]\) есть сетка значений \( x_0 < x_1 < ... < x_n \), в которой известны значения функции \(f\): \(y_0 = f(x_0), y_1 = f(x_1)\), … .

Хотим найти многочлен \(L(x)\) степени не выше \(n\), гарантированно проходящий через точки \(\{ (x_i, y_i) \}_{i=0}^{n}\). Есть надежда, что \(L(x)\) будет близок к \(f(x)\) на всём отрезке \([a, b]\). Назовём \(L(x)\) интерполяционным многочленом. Такая процедура иногда называется Лагранжевой интерполяцией или алгебраической интерполяцией.

Нахождение многочлена \(n\)-й степени, проходящего через \(n+1\) заданную точку - задача со школы (ну или с 1 семестра). Достаточно решить систему линейных уравнений на \(a_i\):

\[\begin{split}\left\{\begin{array}{l} a_{0}+a_{1} x_{0}+a_{2} x_{0}^{2}+\cdots+a_{n} x_{0}^{n}=y_{0} \\ a_{0}+a_{1} x_{1}+a_{2} x_{1}^{2}+\cdots+a_{n} x_{1}^{n}=y_{1} \\ \ldots \ldots \ldots \ldots \ldots \ldots \ldots \ldots \ldots \ldots \ldots \\ a_{0}+a_{1} x_{n}+a_{2} x_{n}^{2}+\cdots+a_{n} x_{n}^{n}=y_{n} \end{array}\right. \end{split}\]

Её можно переписать в матричном виде: \(Xa = y\):

\[\begin{split} \left[\begin{array}{ccccc} 1 & x_{0} & x_{0}^{2} & \ldots & x_{0}^{n} \\ 1 & x_{1} & x_{1}^{2} & \ldots & x_{1}^{n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n} & x_{n}^{2} & \ldots & x_{n}^{n} \end{array}\right]\left[\begin{array}{c} a_{0} \\ a_{1} \\ \vdots \\ a_{n} \end{array}\right]=\left[\begin{array}{c} y_{0} \\ y_{1} \\ \vdots \\ y_{n} \end{array}\right] \end{split}\]

Найдя \(a_i\), итоговый интерполяционный многочлен (aka полиномиальный интерполянт) будет иметь вид:

\[ L_n(x) = a_{0}+a_{1} x +a_{2} x^{2}+\cdots+a_{n} x^{n} \]

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

Играясь с количеством точек и длиной отрезка интерполяции, можно сделать несколько выводов:

• Чем меньше длина отрезка, тем лучше аппроксимация исходной функции и экстраполяция за отрезок

• Для некоторых функций на заданном отрезке увеличение количества точек до определенного предела улучшает аппроксимацию ( \(\sin x,\; \mathcal{N}(x)\)), в то время как для других \(\left(\text{runge}(x)=\cfrac{1}{1+25 \cdot x^2}\right)\) сразу же возникают возрастающие осцилляции.

• Для \(\text{runge}(x)\) при малых \(n\), а для остальных функций при \(n \approx 50\) картина очень напоминает переобучение. И действительно, фактически полиномиальный интерполянт равносилен полиномиальной линейной регрессии с максимальным количеством признаков (гарантированное переобучение). Таким образом, полиномиальная интерполяция автоматически крайне чувствительна к шуму в исходных данных.

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

Самые распространённые формы записи интерполяционного многочлена:

Решать систему выше можно, но абсолютно не нужно. Она решена в общем виде до нас, и интерполяционный многочлен можно написать сразу же, как известны набор точек \(\{ (x_i, y_i) \}_{i=0}^{n}\), правда будет записан он хитро.

Популярны две записи:

• Многочлен Лагранжа — удобен для теоретических исследований, оценок погрешностей, доказательства теорем и аналитических расчетов (в том числе с sympy), но неудобен с точки зрения численных расчётов:

\[\begin{split} \begin{aligned} &L_n(x)=\sum_{i=0}^{n} y_{i} \cdot l_{i}(x) \\ & \\ &l_{i}(x)=\prod_{j=0, \; j \neq i}^{n} \frac{x-x_{j}}{x_{i}-x_{j}}=\frac{x-x_{0}}{x_{i}-x_{0}} \cdots \frac{x-x_{i-1}}{x_{i}-x_{i-1}} \cdot \frac{x-x_{i+1}}{x_{i}-x_{i+1}} \cdots \frac{x-x_{n}}{x_{i}-x_{n}} \end{aligned} \end{split}\]

• Многочлен Ньютона — позволяет легко модифицировать многочлен при добавлении ещё одного узла (в этом случае всего лишь добавляется ещё одно слагаемое).

Процесс записи многочлена Ньютона проходит через разделенную разность — формальное обобщение производной для дискретного набора точек:

\[\begin{split}\begin{aligned} &f\left(x_{0}, \; x_{1}\right)=\frac{f\left(x_{1}\right)-f\left(x_{0}\right)}{x_{1}-x_{0}}, \\ & \\ &f\left(x_{0} ,\; x_{1} ,\; x_{2}\right)=\frac{f\left(x_{1} ,\; x_{2}\right)-f\left(x_{0} ,\; x_{1}\right)}{x_{2}-x_{0}}=\frac{\frac{f\left(x_{2}\right)-f\left(x_{1}\right)}{x_{2}-x_{1}}-\frac{f\left(x_{1}\right)-f\left(x_{0}\right)}{x_{1}-x_{0}}}{x_{2}-x_{0}}, \\ & \\ &f\left(x_{0} ,\; x_{1} ,\; \ldots ,\; x_{n-1} ,\; x_{n}\right)=\frac{f\left(x_{1} ,\; \ldots ,\; x_{n-1} ,\; x_{n}\right)-f\left(x_{0} ,\; x_{1} ,\; \ldots ,\; x_{n-1}\right)}{x_{n}-x_{0}}, \\ & \\ &f\left(x_{j} ,\; x_{j+1} ,\; \ldots ,\; x_{j+k-1} ; x_{j+k}\right)=\frac{f\left(x_{j+1} ,\; \ldots ,\; x_{j+k-1} ,\; x_{j+k}\right)-f\left(x_{j} ,\; x_{j+1} ,\; \ldots ,\; x_{j+k-1}\right)}{x_{j+k}-x_{j}} \end{aligned} \end{split}\]

Заметьте, что при заданном наборе точек \(\{ (x_i, y_i) \}_{i=0}^{n}\) выше написаны какие-то числа, которые можно явно вычислить.

С помощью разделённых разностей можно написать интерполяционный многочлен Ньютона «вперёд» (по нарастанию индексов у точек):

\[ \begin{aligned} L_{n}(x) &=f\left(x_{0}\right)+f\left(x_{0} ,\; x_{1}\right) \cdot\left(x-x_{0}\right)+f\left(x_{0} ,\; x_{1} ,\; x_{2}\right) \cdot\left(x-x_{0}\right) \cdot\left(x-x_{1}\right)+\ldots+f\left(x_{0} ,\; \ldots ,\; x_{n}\right) \cdot \prod_{k=0}^{n-1}\left(x-x_{k}\right) \end{aligned} \]

и аналогично интерполяционный многочлен Ньютона «назад» (по убыванию индексов точек):

\[ \begin{aligned} L_{n}(x) &=f\left(x_{n}\right)+f\left(x_{n} ,\; x_{n-1}\right) \cdot\left(x-x_{n}\right)+f\left(x_{n} ,\; x_{n-1} ,\; x_{n-2}\right) \cdot\left(x-x_{n}\right) \cdot\left(x-x_{n-1}\right)+\ldots+f\left(x_{n} ,\; \ldots ,\; x_{0}\right) \cdot \prod_{k=1}^{n}\left(x-x_{k}\right) \end{aligned} \]

Примечание. С помощью многочлена Ньютона можно также получить следующее представление разделённых разностей в виде отношения определителей:

\[\begin{split}f\left(x_{0} ; \ldots ; x_{n}\right)=\frac{\left|\begin{array}{ccccc} 1 & x_{0} & \ldots & x_{0}^{n-1} & f\left(x_{0}\right) \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 1 & x_{n} & \ldots & x_{n}^{n-1} & f\left(x_{n}\right) \end{array}\right|}{\left|\begin{array}{ccccc} 1 & x_{0} & \ldots & x_{0}^{n-1} & x_{0}^{n} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 1 & x_{n} & \ldots & x_{n}^{n-1} & x_{n}^{n} \end{array}\right|} \end{split}\]

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

Аналитическая интерполяция с помощью sympy

Покажем реализованные python методы, позволяющие нам быстро найти интерполяционный многочлен.

import numpy as np
import matplotlib.pyplot as plt

from matplotlib.pyplot import figure
figure(figsize=(4, 3), dpi=160)

x = np.linspace(0, 1, 10)
y = x**2 - 0.5*x**7

plt.plot(x, y)

[<matplotlib.lines.Line2D at 0x7fa0336dc2d0>]

Имеем 10 точек какого-то многочлена. Построим интерполяционный многочлен 9 степени (Напоминаю, что \(L_n(x)\) строится по \(n+1\) точке) по этому набору точек.

import sympy as smp
from sympy.polys.polyfuncs import interpolate

t = smp.Symbol('t')

points = [(x[i], y[i]) for i in range(len(x))] # Такой формат данных подаётся на вход interpolate

interpolate(points, t)

\[\displaystyle - 6.13908923696727 \cdot 10^{-12} t^{9} + 1.2732925824821 \cdot 10^{-11} t^{8} - 0.500000000152795 t^{7} - 1.16415321826935 \cdot 10^{-10} t^{6} + 7.82165443524718 \cdot 10^{-11} t^{5} - 9.09494701772928 \cdot 10^{-13} t^{4} + 1.13686837721616 \cdot 10^{-12} t^{3} + 1.0000000000004 t^{2} - 7.99360577730113 \cdot 10^{-15} t\]

Видим, что получившийся многочлен очень похож на изначальный - т.е. мы действительно построили интерполяционный многочлен. Мелкие отличия связаны с конечной точностью чисел в x и y.

Построим результат на более мелкой сетке и сравним графики.

f = smp.lambdify(t, interpolate(points, t), 'numpy') # Сделаем python функцию из нашего многочлена

xnew = np.linspace(0, 1, 100)

figure(figsize=(4, 3), dpi=160)

plt.plot(x, y, 'o', xnew, f(xnew), '-')
plt.legend(['data', 'Interpolation'], loc = 'best')
plt.show()

Т.е. мы смогли точно найти все промежуточные значения многочлена по 10 точкам.