TODO

• Убрать вообще пример с матрицей Вандерморда, рисование графиков перенести к fit

Интерполяция

wikipedia

Задача: построить достаточно гладкую кривую, которая гарантированно проходит через заданные точки.

Зачем? Например, проектируем некоторое изделие. Вам нужно задать ограниченное число измерений этого изделия, а по ним нужно составить программу для станка с ЧПУ, регламентирующую положение режущего инструмента с шагом в 10-100-1000 раз меньшим, чем шаг измерений.

Как это решать? Хочется провести кривую из заданного класса, так чтобы она проходила через интересующие нас точки.

Полином какой степени нужен?

Полином степени n - это функция Pn(x) = a_n*xⁿ + a_n-1*x^n-1 + ... + a₂*x² + a₁*x + a₀

Возьмем точку на плоскости. Через нее можно провести единственную прямую вида y = a. Кривых бОльших порядков можно провести неограниченно много (наклонные прямые, параболы и тп).

Возьмем 2 точки на плоскости. Через них можно провести единственную прямую вида y = ax+b. Провести через произвольные 2 точки прямую вида y=a невозможно, а парабол через 2 точки можно провести сколь угодно много.

Для N точек нужно взять интерполяционный полином степени N-1.

Термины

Пусть значения функции φ известны только в этих точках: y_i = φ(x_i), i = 1..N

Задача интерполяции состоит в поиске такой функции F из заданного класса функций , что

F(x_i) = y_i, i=1..N

• точки x_i называют узлы интерполяции
• пары (x_i, y_i) называют точками данных или базовыми точками .
• Разность между соседними значениями dx_i = x_i - x_i-1 это шаг интерполяционной сетки . Он может быть как переменным, так и постоянным.
• Функция F(x) - интерполирующая функция или интерполянт

Базисы функций для интерполяции

Мы будем записывать интерполяционную функцию как линейную комбинацию некоторых базисных функций f_j(x):

F(x) = c₁f₁(x) + c₂f₂(x) +..+ c_Nf_N(x)

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

Import

from scipy import interpolate
from numpy import polynomial as P
import numpy as np
from scipy import linalg
import matplotlib.pyplot as plt

Если мы работаем в notebook и хотим рисовать графики, нужно написать где именно их показывать

%matplotlib inline

Полиномы

Степенные функции 1, x, x², .. x^N - будем пока решать задачу интерполяции на этом базисе.

В numpy есть пакет polynomial, который мы импортировали как Р.

Зададим полином

В примерах задаем полином 1+2x+3x²

Вариант 1: задаем через список коэффициентов а0, а1, ... аn.

In [6]: p1 = P.Polynomial([1, 2, 3])
In [7]: p1
Out[7]: Polynomial([ 1., 2., 3.], [-1, 1], [-1, 1])

Вариант 2: задаем корни уравнения f(x) = 0. Например, x = -1 и x = 1

In [8]: p2 = P.Polynomial.fromroots([-1, 1])
In [9]: p2
Out[9]: Polynomial([-1., 0., 1.], [-1., 1.], [-1., 1.])

Трансформации (не обязательный материал, можно пропустить)

Не обязательный материал. Подробнее смотрите в документации.

Polynomial([-1., 0., 1.], [-1., 1.], [-1., 1.])

В полном описании полинома кроме коэффициентов есть еще два параметра: domain и window. Они задают как с помощью масштабирования и трансформации исходный интервал [domain[0], domain[1]] превращается в результирующий интервал [window[0], window[1]].

Коэффициенты, исходный и результирующий интервал можно получить из полинома с помощью методов:

In [12]: p1.coef
Out[12]: array([ 1., 2., 3.])
In [13]: p1.domain
Out[13]: array([-1, 1])
In [14]: p1.window
Out[14]: array([-1, 1])

Найдем корни полинома roots()

In [10]: p1.roots()
Out[10]: array([-0.33333333-0.47140452j, -0.33333333+0.47140452j])
In [11]: p2.roots()
Out[11]: array([-1., 1.])

Вычислим значения полинома в точках

Значение полинома в 1 точке:

In [14]: p1(1.5)
Out[14]: 10.75

Найдем значения полинома в нескольких точках. Например, в точках 1.5, 2.5, 3.5.

Просто вызовем p1 от списка точек и получим массив значений.

In [15]: p1([1.5, 2.5, 3.5])
Out[15]: array([ 10.75, 24.75, 44.75])

Арифметические операции над полиномами

Для полиномов определены арифметические операции +, -, *, / и так далее.

Операция // используется для полиномиального деления.

Пусть один полином будет p1(x) = (x-3)(x-2)(x-1), другой p2(x) = (x-2). Тогда деление одного полинома на другой будет p3 = (x-3)(x-2)(x-1)/(x-2) = (x-3)(x-1)

Убедимся, что это так. Зададим полиномы через их корни. Разделим один на другой и найдем корни полученного полинома.

In [16]: p1 = P.Polynomial.fromroots([1, 2, 3])
In [17]: p1
Out[17]: Polynomial([ -6., 11., -6., 1.], [-1., 1.], [-1., 1.])
In [18]: p2 = P.Polynomial.fromroots([2])
In [19]: p2
Out[19]: Polynomial([-2., 1.], [-1., 1.], [-1., 1.])
In [20]: p3 = p1 // p2
In [21]: p3
Out[21]: Polynomial([ 3., -4., 1.], [-1., 1.], [-1., 1.])
In [22]: p3.roots()
Out[22]: array([ 1., 3.])

Производная

polyder(p, m=1) взять производную полинома степени m (m>0)

Пример из документации:

>>> from numpy.polynomial import polynomial as Pl   # тут другой import
>>> c = (1,2,3,4)             # 1 + 2x + 3x**2 + 4x**3
>>> Pl.polyder(c)             # (d/dx)(c) = 2 + 6x + 12x**2
array([  2.,   6.,  12.])
>>> Pl.polyder(c,3)           # (d**3/dx**3)(c) = 24
array([ 24.])
>>> Pl.polyder(c,scl=-1)      # (d/d(-x))(c) = -2 - 6x - 12x**2
array([ -2.,  -6., -12.])
>>> Pl.polyder(c,2,-1)        # (d**2/d(-x)**2)(c) = 6 + 24x
array([  6.,  24.])

Заметим, что polyder принимает аргументом не полином, а список коэффициентов. То есть для полинома нужно писать так:

In [23]: p = P.Polynomial([1,2,3])
In [24]: P.polynomial.polyder(p.coef)  # возьмем из полинома список коэффициентов
Out[24]: array([ 2.,  6.])             # d(1+2x+3x**2)/dx = 2+6x

Интеграл

m=1, k=[], lbnd=0, scl=1, axis=0) взять m раз первообразную полинома с коэффициентами с, добавляя каждый раз константу из списка k и умножая на scl.

Для интегрирования 1 + 2*x + 3*x**2 возьмите c = [1,2,3].

Для интегрирования 1 + 2*x + 3*y + 4*x*y возьмите c = 1,2?, для интегрирования по x возьмите axis=0, а для интегрирования по y укажите axis=1.

На каждой итерации полученная первообразная умножается на scl и к ней прибавляется константа k.

Если m≥1, то умножение на scl делается m раз. То есть, если нужно сделать линейное преобразование u = ax+v, то du = a*dx, и dx = du/a. Это 1/a нужно указать в scl=1/a.
k тоже можно указать списком размера m или меньше. Тогда константы будут добавляться на каждом этапе интегрирования.

Пример из документации:

>>> from numpy.polynomial import polynomial as P
>>> c = (1,2,3)
>>> P.polyint(c)                  # Integral(1+2x+3x**2)dx = x + x**2 + x**3, то есть array([0, 1, 1, 1])
array([ 0.,  1.,  1.,  1.])
>>> P.polyint(c,3)                # 3 раза взяли интеграл, должны получить array([0, 0, 0, 1/6, 1/12, 1/20])
array([ 0. ,  0. ,  0. ,  0.16666667,  0.08333333, 0.05 ])
>>> P.polyint(c,k=3)              # к первообразной прибавили k=3, получили array([3, 1, 1, 1])
array([ 3.,  1.,  1.,  1.])
>>> P.polyint(c,lbnd=-2)          # к первообразной прибавили ее значение при x=-2, получили array([6, 1, 1, 1])
array([ 6.,  1.,  1.,  1.])
>>> P.polyint(c,scl=-2)           # первообразную умножили на scl=-2, получили array([0, -2, -2, -2])
array([ 0., -2., -2., -2.])

Заметим, что polyint принимает аргументом не полином, а список коэффициентов. То есть для полинома нужно писать так:

In [23]: p = P.Polynomial([1,2,3])
In [24]: P.polynomial.polyint(p.coef)  # возьмем из полинома список коэффициентов
Out[24]: array([0., 1., 1., 1.])       # Integral(1+2x+3x**2)dx = x + x**2 + x**3 + C, где С взяли равным 0.

Другие типы многочленов

Аналогично можем задать многочлен Чебышева с коэффициентами 1, 2, 3, т.е. полином 1T₁+2T₂+3T₃, где T_i(x) - многочлен Чебышева порядка i

In [23]: c1 = P.Chebyshev([1, 2, 3])
In [24]: c1
Out[24]: Chebyshev([ 1., 2., 3.], [-1, 1], [-1, 1])

вычислить его корни

In [25]: c1.roots()
Out[25]: array([-0.76759188, 0.43425855])

Все многочлены других классов имеют те же методы, что и класс Polynomial, который описан выше. Используют их тоже аналогично.

Задать многочлены из корней:

In [26]: c1 = P.Chebyshev.fromroots([-1, 1])
In [27]: c1
Out[27]: Chebyshev([-0.5, 0. , 0.5], [-1., 1.], [-1., 1.])
In [28]: l1 = P.Legendre.fromroots([-1, 1])
In [29]: l1
Out[29]: Legendre([-0.66666667, 0. , 0.66666667], [-1., 1.], [-1., 1.])

Подсчитать значения многочленов в точках:

In [30]: c1([0.5, 1.5, 2.5])
Out[30]: array([-0.75, 1.25, 5.25])
In [31]: l1([0.5, 1.5, 2.5])
Out[31]: array([-0.75, 1.25, 5.25])

Интерполяция алгебраическими многочленами

Рассмотрим набор точек с координатами (1, 1), (2, 3), (3, 5) и (4, 4). Создадим соответствующие NumPy? массивы для х и у координат:

In [32]: x = np.array([1, 2, 3, 4])
In [33]: y = np.array([1, 3, 5, 4])

Для интерполяции многочленом через эти точки, надо использовать полином 3 степени (количество точек - 1). Для интерполяции в степенном базисе найдем полином в этом базисе.

Максимальная степень полинома:

In [34]: deg = len(x) - 1

Метод fit

Используем метод fit для степенного полинома и полинома Чебышева:

In [35]: f1b = P.Polynomial.fit(x, y, deg)
In [36]: f1b
Out[36]: Polynomial([ 4.1875, 3.1875, -1.6875, -1.6875], [ 1., 4.], [-1., 1.])

In [37]: f2b = P.Chebyshev.fit(x, y, deg)
In [38]: f2b
Out[38]: Chebyshev([ 3.34375 , 1.921875, -0.84375 , -0.421875], [ 1., 4.], [-1., 1.])

Получен вектор коэффициентов [ 4.1875, 3.1875, -1.6875, -1.6875] в степенном базисе, то есть интерполяционный полином f (x) = 4.1875 + 3.1875*x - 1.6875 * x² - 1.687 5* x³.

Используем найденный массив коэффициентов для дальнейшего вычисления значения в точке:

In [39]: f1b(2.5)
Out[39]: 4.1875000000000018

Как и ожидалось, в другом базисе получились другие коэффициенты и интерполяционный многочлен в базисе Чебышева f (x) = 3.34375*T₀(x)- 1.921875*T₁(x)- 0.84375*T₂(x) - 0.421875*T₃(x) Однако, не зависимо от базиса, интерполяционный многочлен единственный и вычисляя значение в точке получаем такое же число:

In [40]: f2b(2.5)
Out[40]: 4.1875000000000018

Посмотрим на эти два полинома вместе на одном графике:

In [45]: xx = np.linspace(x.min(), x.max(), 100) # supersampled [x[0], x[-1]] interval
In [45]: fig, ax = plt.subplots(1, 1, figsize=(12, 4))
...: ax.plot(xx, f1b(xx), 'b', lw=2, label='Power basis interp.')
...: ax.plot(xx, f2b(xx), 'r--', lw=2, label='Chebyshev basis interp.')
...: ax.scatter(x, y, label='data points')
...: ax.legend(loc=4)
...: ax.set_xticks(x)
...: ax.set_ylabel(r"$y$", fontsize=18)
...: ax.set_xlabel(r"$x$", fontsize=18)

Если вы пишите код в отдельном файле, а не в notebook, не забудьте дописать, чтобы показать график:

plt.show()

fit: domain и window, погрешность интерполяции

Заметим, что этот метод автоматически устанавливает атрибут domain в правильное значение (в примере отрезок был от 1 до 4), коэффициенты устанавливаются тоже соответствующим образом.

Отображение интерполяционных данных в отрезок, который лучше подходит для выбранного базиса, может существенно улучшить численную стабильность интерполяции. Например, использование полиномов Чебышева на отрезке [-1,1] а не на заданном [1,4] уменьшает число обусловленности с почти 4660 до 1.85:

Чем ближе число обусловленности к 1, тем меньше погрешность интерполянта.

In [50]: np.linalg.cond(P.chebyshev.chebvander(x, deg))     # x = [1, 2, 3, 4]
Out[50]: 4659.7384241399586
In [51]: np.linalg.cond(P.chebyshev.chebvander((2*x-5)/3.0, deg))  # (2*x-5)/3.0 = [-1, -5/3, 1/3, 1]
Out[51]: 1.8542033440472896

откуда взялись эти х и (2*x-5)/3.0? Уравнение трансформации, переводящее 1 в -1 и 4 в 1 есть (2*x-5)/3.0 (его можно найти самим)

Больше степень - больше биения. Пример интерполяции функции Рунге

Интерполяция многочленами небольшого количества точек - мощный и полезный математический инструмент. Но если количество точек растет, то растет и наибольшая степень интерполяционного многочлена, что приводит к проблемам.

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

Возьмем функцию Рунге f (x) = 1/(1+25x^2), достаточно ровную на интервале [-1,1]. Результат интерполяции почти расходится в промежутке между точками у конца интервала интерполяции.

Для иллюстрации этого поведения напишем функцию runge, которая реализует функцию Рунге, и функцию runge_interpolate, которая находит интерполяцию n-той степени в степенном базисе.

In [52]: def runge(x):
...: return 1/(1 + 25 * x**2)
In [53]: def runge_interpolate(n):
...: x = np.linspace(-1, 1, n)
...: p = P.Polynomial.fit(x, runge(x), deg=n)
...: return x, p

Нарисуем графики исходной функции и интерполяционного многочлена 13 и 14 порядков:

In [54]: xx = np.linspace(-1, 1, 250)
In [55]: fig, ax = plt.subplots(1, 1, figsize=(8, 4))
...: ax.plot(xx, runge(xx), 'k', lw=2, label="Runge function")
...: # 13th order interpolation of the Runge function
...: n = 13
...: x, p = runge_interpolate(n)
...: ax.plot(x, runge(x), 'ro')
...: ax.plot(xx, p(xx), 'r', label='interp. order %d' % n)
...: # 14th order interpolation of the Runge function
...: n = 14
...: x, p = runge_interpolate(n)
...: ax.plot(x, runge(x), 'go')
...: ax.plot(xx, p(xx), 'g', label='interp. order %d' % n)
...: ax.legend(loc=8)
...: ax.set_xlim(-1.1, 1.1)
...: ax.set_ylim(-1, 2)
...: ax.set_xticks([-1, -0.5, 0, 0.5, 1])
...: ax.set_ylabel(r"$y$", fontsize=18)
...: ax.set_xlabel(r"$x$", fontsize=18)

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

Интерполяция сплайнами

Для набора из n точек {x_i, y_i}, на интервале [x₀, x_n-1] получим n-1 подинтервал [x_i, x_i+1]. Внутренняя точка, которая соединяет два подинтервала, будем называть узлом в терминах кусочно-полиномиальной интерполяции. Для интерполяции n точек, будем использовать кусочное полиномы степени k на каждом подинтервале, т.е. мы должны определить (k+1)(n-1) неизвестных параметра. Значения в узлах дают 2(n-1) уравнения. Этих уравнений недостаточно для нахождения всех неизвестных. Но требование непрерывности производных более высоких порядков дадут нам необходимые уравнения.

Это условие гарантирует, что кусочно-полиномиальный интерполянт будет достаточно гладким.

Сплайн - это специальный тип кусочно-полиномиального интерполянта: кусочный полином степени k является сплайном, если он непрерывно дифференцируем k-1 раз.

Наиболее популярным является сплайн 3 степени, k=3, который требует 4(n-1) параметр. В этом случае непрерывность двух производных в n-2 точках дает 2(n-2) дополнительных уравнения, повышая общее количество уравнений до 2 (n -1)+ 2(n - 2) = 4(n -1)- 2.

Остается 2 неопределенных параметра, которые мы должны определить другим способом.

Обычно задают дополнительное требование, чтобы производная второго порядка на концах точек была равна 0 (превращая в настоящий сплайн). Это дает еще 2 уравнения.

Модуль интерполяции SciPy? предоставляет несколько функций и классов для представления сплайн-интерполяции. Например, мы можем использвать функцию interpolate.interp1d, которая берет массивы х и у для точек как первый и второй аргументы. Дополнительные аргументы позволяют задать тип и порядок интерполяции. Обычно задают kind=3 (или, что тоже самое, kind='cubic') для вычисления кубической интерполяции. Эта функция возвращает класс объектов, которые можно вызывать как функцию и которая может вычислять значения для различных значений х, используя вызовы функции.

Альтернативная сплайн-функция - это interpolate.InterpolatedUnivariateSpline, которая так же принимает массивы х и у в качестве первых двух аргументов, но которая использует аргумент k (вместо kind) для задания порядка сплайнов.

Посмотрим, как можно использовать функцию interpolate.interp1d для фунцкии Рунгде и возьмем сплайны 3 степени. Для этого мы сначала создадим массивы NumPy? для х и у координат в тех же точках. Дальше вызовем interpolate.interp1d с k=3 для получения кубических сплайнов на заданных данных.

In [56]: x = np.linspace(-1, 1, 11)
In [57]: y = runge(x)
In [58]: f_i = interpolate.interp1d(x, y, kind=3)

Для того, чтобы посмотреть насколько хороши сплайны, нарисуем исходную функцию Рунге, точки и кубические сплайны.

In [59]: xx = np.linspace(-1, 1, 100)
In [60]: fig, ax = plt.subplots(figsize=(8, 4))
...: ax.plot(xx, runge(xx), 'k', lw=1, label="Runges function")
...: ax.plot(x, y, 'ro', label='sample points')
...: ax.plot(xx, f_i(xx), 'r--', lw=2, label='spline order 3')
...: ax.legend()
...: ax.set_xticks([-1, -0.5, 0, 0.5, 1])
...: ax.set_ylabel(r"$y$", fontsize=18)
...: ax.set_xlabel(r"$x$", fontsize=18)

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

Для иллюстрации влияния порядка сплайна, решим задачу интерполяции на точках (0,3), (1, 4), (2, 3.5), (4, 2), (5, 1.5), (6, 1.25) и (7, 0.7) со сплайнами возрастающей степени.

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

In [61]: x= np.array([0, 1, 2, 3, 4, 5, 6, 7])
In [62]: y= np.array([3, 4, 3.5, 2, 1, 1.5, 1.25, 0.9])
In [63]: xx = np.linspace(x.min(), x.max(), 100)
In [64]: fig, ax = plt.subplots(figsize=(8, 4))
...: ax.scatter(x, y)
...:
...: for n in [1, 2, 3, 6]:
...:     f = interpolate.interp1d(x, y, kind=n)
...:     ax.plot(xx, f(xx), label='order %d' % n)
...:
...: ax.legend()
...: ax.set_ylabel(r"$y$", fontsize=18)
...: ax.set_xlabel(r"$x$", fontsize=18)

Можно проще?

Можно. Если точек немного используйте интерполяцию (степенную) модуля sympy:

http://docs.sympy.org/latest/modules/polys/reference.html

from sympy.polys.polyfuncs import interpolate
import sympy
x = sympy.symbols('x')

Для точек {(1,1), (2, 4), (3, 9), (4, 16)}

>>> interpolate([1, 4, 9, 16], x)
x**2

или явно задайте список координат:

>>> interpolate([(1, 1), (2, 4), (3, 9)], x)
x**2

х и у координаты могут задаваться как ключ и значение словаря.

>>> interpolate([(-1, 2), (1, 2), (2, 5)], x)
x**2 + 1
>>> interpolate({-1: 2, 1: 2, 2: 5}, x)
x**2 + 1

Для конструирования полиномов Лагранжа по n точкам используйте

sympy.polys.specialpolys.interpolating_poly(n, x, X='x', Y='y')[source]

Задание

Есть набор точек:

X:[1,2,3,4,5,6,7];

Y:[2,0,-1,1,2,3,2];

Построить по этим точкам

• Степенной многочлен
• Многочлен в форме Лагранжа (или Чебышева)
• Кусочно-линейный интерполянт степени 1
• Кусочно-линейный интерполянт степени 3
• Отобразить на одном графике.

Дополнительные вопросы

* Для получившихся интерполянтов найти максимум на отрезке от 1 до 7.

* Для получившихся интерполянтов найти максимальное по всем парам соседних точек отклонение значения интерполянта от значения в точке пары с наименьшей x-координатой.

Это читать уже не нужно, старая теория (кому интересно)

Интерполяция алгебраическими многочленами

Теория

Через вычисление матрицы Вандермонда - НЕ НАДО ТАК РЕШАТЬ! ИСПОЛЬЗУЙТЕ fit, ОН ПРОЩЕ

Классы многочленов, описанные выше, дают полезные функции для интерполяции многочленами. Например, решим линейное уравнение для задачи интерполяции F(x)c = y, где х и у - вектора с координатами xi и yi. Для решения задачи интерполяции мы должны вычислить матрицу F(x) в выбранном базисе и далее решить систему линейных уравнений. Каждый класс многочленов имеет функции для вычисления (обобщенной) матрицы Вандермонда для соответствующего базиса.

Например, для степенного базиса мы используем np.polynomial.polynomial.polyvander, для базиса из многочленов Чебышева берем функцию np.polynomial.chebyshev.chebvander и так далее.

Полный список функций приведен в документации по модулю np.polynomial и его подмодулям.

Используя эти функции для создания матриц Вандерморта, мы можем легко представить интерполяционный многочлен в любом базисе.

Например, рассмотрим набор точек с координатами (1, 1), (2, 3), (3, 5) и (4, 4). Создадим соответствующие NumPy? массивы для х и у координат:

In [32]: x = np.array([1, 2, 3, 4])
In [33]: y = np.array([1, 3, 5, 4])

Для интерполяции многочленом через эти точки, надо использовать полином 3 степени (количество точек - 1). Для интерполяции в степенном базисе найдем коэффициенты c_i такие что f(x) = c₁x⁰+c₂x¹+c₃x²+c₄x³ Чтобы найти эти коэффициенты, вычислим матрицу Вандермонда и решим систему интерполяционных уравнений:

In [34]: deg = len(x) - 1
In [35]: A = P.polynomial.polyvander(x, deg)
In [36]: c = linalg.solve(A, y)
In [37]: c
Out[37]: array([ 2. , -3.5, 3. , -0.5])

Получен вектор коэффициентов [2, -3.5, 3, -0.5], то есть интерполяционный полином f (x) = 2 - 3.5x + 3x² - 0.5x³.

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

In [38]: f1 = P.Polynomial(c)
In [39]: f1(2.5)
Out[39]: 4.1875

Найдем интерполяционный многочлен в другом базисе, например, многочленов Чебышева:

In [40]: A = P.chebyshev.chebvander(x, deg)
In [41]: c = linalg.solve(A, y)
In [42]: c
Out[42]: array([ 3.5 , -3.875, 1.5 , -0.125])

Как и ожидалось, в другом базисе получились другие коэффициенты и интерполяционный многочлен в базисе Чебышева f (x) = 3.5T₀(x)- 3.875T₁(x)+1.5T₂(x)- 0.125T₃(x) Однако, не зависимо от базиса, интерполяционный многочлен единственный и вычисляя значение в точке получаем такое же число:

In [43]: f2 = P.Chebyshev(c)
In [44]: f2(2.5)
Out[44]: 4.1875

Посмотрим на эти два полинома вместе на одном графике:

In [45]: xx = np.linspace(x.min(), x.max(), 100) # supersampled [x[0], x[-1]] interval
In [45]: fig, ax = plt.subplots(1, 1, figsize=(12, 4))
...: ax.plot(xx, f1(xx), 'b', lw=2, label='Power basis interp.')
...: ax.plot(xx, f2(xx), 'r--', lw=2, label='Chebyshev basis interp.')
...: ax.scatter(x, y, label='data points')
...: ax.legend(loc=4)
...: ax.set_xticks(x)
...: ax.set_ylabel(r"$y$", fontsize=18)
...: ax.set_xlabel(r"$x$", fontsize=18)

Это не очень удобно. Хотим увеличить порядок полинома - нужно переписывать уравнение или программировать на эту тему какую-то автоматику. Может, есть что готовое?

Уроки прошлых лет в виде тетради notebook

• interpolate.ipynb: урок по интерполяции 2019 года

Attachment	Action	Size	Date	Who	Comment
fig1.png	manage	23.8 K	22 Apr 2016 - 22:45	TatyanaDerbysheva	рис. 7-1
fig2.png	manage	47.6 K	22 Apr 2016 - 23:31	TatyanaDerbysheva	рис. 7-2
fig3.png	manage	30.2 K	23 Apr 2016 - 00:18	TatyanaDerbysheva	рис. 7-3
fig4.png	manage	35.6 K	23 Apr 2016 - 00:18	TatyanaDerbysheva	рис. 7-4
interpolate.ipynb	manage	164.2 K	27 Nov 2019 - 09:05	TatyanaDerbysheva	урок по интерполяции 2019 года