In [1]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate
Допустим мы знаем значения $f_k=f(x_k)$ некоторой функции $f(x)$ только на некотором множестве аргументов $x_k\in\mathbb R$, $k=1..K$. Мы хотим вычислять $f$ в точках $x$ лежащих между узлами интерполяции $x_k$, для чего мы строим интерполирующую функцию $p(x)$, определенную при всех значениях $x$ и совпадающую с $f$ в узлах интерполяции $f(x)=p(x)$, $k=1..K$. Распространенным (но не единственным) выбором для функции $p$ является многочлен степени $K-1$:
$$p(x)=\sum_{n=0}^{K-1} p_n x^n,$$число коэффициентов которого равно число известных значений функции, что позволяет однозначно найти эти коэффициенты. Формально коэффициенты находятся из системы уравнений:
$$p(x_k)=f_k=\sum_{n=0}^{K-1} p_n x_k^n,$$или в матричном виде $MP=F$, где
$$ F=\begin{pmatrix}f_1\\\vdots\\f_K\end{pmatrix},\quad P=\begin{pmatrix}p_0\\\vdots\\p_{K-1}\end{pmatrix},\quad M=\begin{pmatrix} x_1^0 & \cdots & x_1^{K-1} \\ \vdots & \ddots & \vdots \\ x_K^0 & \cdots & x_K^{K-1} \\ \end{pmatrix}. $$Матрица $M$ называется матрицей Вандермонда. Попробуем построить интерполяционный многочлен таким способом.
In [2]:
# Numpy уже имеет класс для полиномов: numpy.poly1d.
# Для полноты изложения мы реализуем свой класс.
class Poly():
def __init__(self, pn):
"""
Создает многочлен с коэффициентами pn:
self(x) = sum_n pn[n] * x**n.
Коэффициенты pn перечисляются в порядке возрастания степени одночлена.
"""
self.pn = pn
def __call__(self, x):
"""
Вычисляет многочлен на векторе значений x.
"""
a = 1. # Здесь мы накапливаем степени x**n.
p = 0. # Сюда мы помещаем сумму одночленов.
for pn in self.pn:
p += a*pn # Учитываем очередной одночлен.
a *= x # Повышаем степень одночлена
return p
# Вспомогательная функция для счета матрицы Вандермонда.
def vandermonde(xn):
return np.power(xn[:,None], np.arange(len(xn))[None,:])
# Напишем функцию, которая будет находить интерполяционных многочлен через решение системы.
def interp_naive(xn, fn):
"""
Возвращает интерполяционный многочлен, принимающий в точках xp значение fp.
"""
M = vandermonde(xn)
# Мы используем функцию numpy для решения линейных систем.
# Методы решения линейных систем обсуждаются в другой лабораторной работе.
pn = np.linalg.solve(M, fn)
return Poly(pn)
# Возьмем логарифмическую решетку на интервале [1E-6,1] и посмотрим, насколько точно мы восстанавливаем многочлен.
N = 8 # Число узлов интерполяции.
x = np.logspace(-6,0,N) # Точки равномерной решетки
# Будет интерполировать многочлен степени N-1, который зададим случайными коэффициентами.
f = Poly(np.random.randn(N))
y = f(x) # Значения многочлена на решетке.
p = interp_naive(x, y) # Строим интерполяционный многочлен.
z = p(x) # Значения интерполяционного многочлена на решетке.
print("Absolute error of values", np.linalg.norm(z-y))
print("Absolute error of coefficients", np.linalg.norm(f.pn-p.pn))
Построенный интерполяционный многочлен принимает близкие значения к заданным в узлах значениям. Но хотя значения в узлах должны задавать многочлен однозначно, интерполяционный и интерполируемый многочлен имеют значительно отличающиеся коэффициенты. Значения многочленов между узлами также значительно отличаются.
In [3]:
t = np.linspace(0,1,100)
_, ax = plt.subplots(figsize=(10,7))
ax.plot(t, f(t), '-k')
ax.plot(t, p(t), '-r')
ax.plot(x, f(x), '.')
ax.set_xlabel("Argument")
ax.set_ylabel("Function")
ax.legend(["$f(x)$", "$p(x)$"])
plt.show()
__call__, так чтобы он реализовывал схему Горнера. Чем эта схема лучше?На практике интерполяционный многочлен обычно находится в форме многочлена Лагранжа:
$$ p(x)=\sum_{k=1}^K f_k L_k(x),\;\text{где}\; L_k(x)=\prod_{j\neq k}\frac{x-x_j}{x_k-x_j}. $$Для ускорения вычисления многочлена Лагранжа используется схема Эйткена, основанная на рекурсии. Обозначим $p_{i,\ldots,j}$ многочлен Лагранжа, построенный по узлам интерполяции $(x_i,f_i),\ldots,(x_j,f_j)$, в частности искомое $p=p_{1,\ldots,K}$. Справедливо следующее соотношение, выражающие интеполяционный многочлен через такой же с меньшим числом узлов:
$$ p_{i,\ldots,j}=\frac{(x-x_i)p_{i+1,\ldots,j}-(x-x_j)p_{i,\ldots,j-1}(x)}{x_j-x_i}. $$База рекурсии задается очевидным равенством $p_{i}(x)=f_i$.
Интерполяционные формулы Ньютона дают другой популярный способ записи интерполяционного многочлена:
$$ p(x)=\sum_{k=1}^K[x_1,\ldots,x_k]f\prod_{i=1}^k(x-x_j), $$где разделенные разности $[\ldots]f$ заданы рекурсивно:
$$ [x_1,\ldots,x_k,x]f=\frac{[x_1,\ldots,x_{k-1},x]f-[x_1,\ldots,x_{k-1},x_k]f}{x-x_k}. $$Данные формулы можно рассматривать как дискретный вариант формулы Тейлора. На основе формул Ньютона разработан алгоритм Невилла для вычисления интерполяционного многочлена, по существу эквивалентный схеме Эйткена.
scipy.interpolate.lagrange. В документации отмечается, что метод численно неустойчив. Что это означает?Рассмотрим теперь насколько хорошо интерполяционный многочлен прилижает интерполируемую функцию.
In [4]:
# В качестве интерполируемой функции возьмем f(x)=x sin(2x.
def f(x):
return x*np.sin(2*x)
# Будем интерполировать функцию на интервале [x0-r,x0+r], где
x0 = 10
r = 1
# В качестве узлов интерполяции возьмем равномерную решетку из N узлов.
N = 5
xn = np.linspace(x0-r, x0+r, N)
# Построим интерполяционный многочлен.
p = interp_naive(xn, f(xn))
# Оценим точность приближения функции многочленом как максимум
# отклонения значений многочлена от значений функции на интервале.
# Так как мы не можем рассмотреть все точки, то ограничимся
# плотной решеткой.
tn = np.linspace(x0-r, x0+r, 10000)
error = np.abs(f(tn)-p(tn))
print("Error", np.max(error))
_, ax = plt.subplots(figsize=(10,5))
ax.semilogy(tn, error)
ax.set_xlabel("Argument")
ax.set_ylabel("Absolute error")
plt.show()
Найдите погрешность прилижения функции $f$ интерполяционным многолченом $p$ для $x0=10, 100, 1000$ и для $N=5, 10, 15$. Объясните получающиеся результаты.
Постройте график зависимости ошибки от числа узлов интерполяции $N$ для $x0=100$ и $r=5$ в диапазоне $5\leq N \leq 50$.
Повторите задания 9 и 10 для узлов интерполяции Чебышева:
Сравните распределение ошибки внутри интервала $x\in[x0-r,x0+r]$ для равномерно расположенных узлов и для узлов Чебышева.
Повторите задания 9 и 10 для функции $f(x)=|x-1|$, $x0=1$, $r=1$. Объясните наблюдающиеся различия.
Использование интерполяционного полинома очень высокой степени часто приводит к тому, что в некоторых точках погрешность приближения оказывается очень большой. Вместо одного многочлена высокой степени, приближающего функцию на всем интервале, можно использовать несколько многочленов меньше степени, каждый из которых приближает функцию только на подинтервале. Если функция обладает некоторой степенью гладкости, например, несколько ее производных непрерывные функции, то такую же гладкость естественно требовать от результирующего семейства интерполяционных многочленов, что накладывает ограничения на их коэффициенты. Получающаяся кусочно-полиномиальная функция называется сплайном. Кубическим сплайном дефекта 1 называется функция, которая:
Для функции из задания 9 постройте кубический сплайн дефекта 1 с узлами из задания 9. Можете воспользоваться функциями scipy.interpolate.splrep и scipy.interpolate.splev или реализовать свои аналоги.
Изучите зависимость погрешности приближения функции сплайном от числа узлов интерполяции. Сравните с результатом из задания 10. Когда погрешности совпадут?
Как можно обобщить изученные методы интерполяции на кривые в многомерном пространстве?
Как можно интерполировать функции нескольких переменных?
Какие еще способы интерполяции существуют?
In [ ]: