# 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))

"""
Explanation: Интерполяция многочленами
Допустим мы знаем значения $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$ называется матрицей Вандермонда.
Попробуем построить интерполяционный многочлен таким способом.
End of explanation
"""


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()

"""
Explanation: Построенный интерполяционный многочлен принимает близкие значения к заданным в узлах значениям.
Но хотя значения в узлах должны задавать многочлен однозначно, интерполяционный и интерполируемый многочлен имеют значительно отличающиеся коэффициенты. Значения многочленов между узлами также значительно отличаются.
End of explanation
"""


# В качестве интерполируемой функции возьмем 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()

"""
Explanation: Задания

Измените метод __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}.
$$
Данные формулы можно рассматривать как дискретный вариант формулы Тейлора. 
На основе формул Ньютона разработан алгоритм Невилла для вычисления интерполяционного многочлена, по существу эквивалентный схеме Эйткена.
Задания

Реализуйте метод Эйткена вычисления интерполяционного многочлена.  
Если мы попытаемся восстановить многочлен через его значения в точках, аналогично заданию 2, получим ли мы с помощью метода Эйткена ответ точнее, чем через решение системы?
Scipy содержит готовую реализацию интерполяционного многочлена Лагранжа scipy.interpolate.lagrange. В документации отмечается, что метод численно неустойчив. Что это означает?
Ошибки в исходных данных для построения интерполяционного многочлена вызывают ошибки при вычислении интерполяционного многочлена в промежуточных точках. При каком расположении узлов интерполяция многочленом Лагранжа имеет наименьшую ошибку? Как это связано с численной устойчивостью? 

Рассмотрим теперь насколько хорошо интерполяционный многочлен  прилижает интерполируемую функцию.
End of explanation
"""


