<jupyter_start><jupyter_text>Занятие 5 Прикладная алгебра и численные методы Интерполяция: полиномы Лагранжа, сплайны и кривые Безье<jupyter_code>import numpy as np
import scipy.linalg
import sympy
from sympy import S
import matplotlib.pyplot as plt
from copy import deepcopy
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Многочлен Лагранжа$$f(x) = L(x) = \sum_{i=1}^n y_i \frac{(x - x_0)(x - x_1)...(x - x_{i-1})(x - x_{i+1})...(x - x_n)}{(x_i - x_0)(x_i - x_1)...(x_i - x_{i-1})(x_i - x_{i+1})...(x_i - x_n)}, \quad y_i = f(x_i), \quad i = 0, ..., n.$$ Пример 1Построить параболу по трем точкам (-1, -2), (0, -1), (1, 2).<jupyter_code>x = S('x')
X = (-1, 0, 1)
Y = (-2, -1, 2)
L = 0
for i in range(3):
    Li = Y[i]
    for j in range(3):
        if i != j:
            Li *= (x - X[j])/(X[i] - X[j])
    L += Li
display(L, sympy.simplify(sympy.expand(L)))<jupyter_output><empty_output><jupyter_text>Пример 2 Линейная регрессияДля данных Примера 1 построим параболу методом линейной регрессии (исключительно для тренировки и повторения, в жизни так не делают!)Составим матрицу  из $A$ из значений $x$ и $x^2$ и матрицу $Y$, далее все как обычно - псевдообратная матрица и произведение псевдообратной к  $A$ и  $Y$:<jupyter_code>X2 = [ -1, 0, 1]
A2 = np.array([[1, item, item**2] for item in X2])
Y2 = np.array([[-2], [-1], [2]])
res2 = np.matmul(np.linalg.pinv(A2), Y2)
a2, b2, c2 = [round(item, 2) for item in  res2[:, 0]]
display('A', A2,'Y', Y)
print('a =', a2, 'b =', b2, 'c =', c2)<jupyter_output><empty_output><jupyter_text>Изобразим на графике точки исходных данных и полученную линию регрессии:<jupyter_code>X_ls2 = np.linspace(X2[0], X2[-1])
plt.plot(X_ls2, a2 + b2*X_ls2 + c2*X_ls2**2, label='$y = {a} + {b}x + {c}x^2$'.format(a=round(a2), 
                                                                                      b=round(b2), c=round(c2)))
plt.scatter(X2, Y2, color='red', label='data')
plt.grid()
plt.legend()<jupyter_output><empty_output><jupyter_text>Пример 3 Линейная регрессияДля данных точек построим наилучшую в смысле среднеквадратического отклонения параболу методом линейной регрессии.$$\begin{matrix}X & -3 & -1 & 0 & 1 & 3\\Y & -4 & -0.8 & 1.6 & 2.3 & 1.5\end{matrix}$$<jupyter_code>X3 = [ -3, -1, 0, 1, 3]
A3 = np.array([[1, item, item**2] for item in X3])
Y3 = np.array([[-4], [-0.8], [1.6], [2.3], [1.5]])
res3 = np.matmul(np.linalg.pinv(A3), Y3)
a3, b3, c3 = [round(item, 2) for item in  res3[:, 0]]
print('a =', a3, 'b =', b3, 'c =', c3)<jupyter_output>a = 1.23 b = 0.98 c = -0.28<jupyter_text>Построим график, отметим на нем исходные точки.<jupyter_code>X_ls = np.linspace(X3[0], X3[-1])
plt.plot(X_ls, a3 + b3*X_ls + c3*X_ls**2)
plt.scatter(X3, Y3, color='red')
plt.grid()<jupyter_output><empty_output><jupyter_text>Пример 4Построим полином Лагранжа для данных Примера 3:<jupyter_code>x = S('x')
Y3 = [-4, -S(8)/10, S(16)/10, S(23)/10, S(15)/10]
L3 = 0
n = len(X3)
for i, Li in enumerate(Y3):
    for j in range(n):
        if i != j:
            Li *= (x - X3[j])/(X3[i] - X3[j])
    L3 += Li
L3 = sympy.simplify(sympy.expand(L3))    
display(L3)
X_ls = np.linspace(X3[0], X3[-1])
Y_regr = [L3.subs(x, item) for item in X_ls]
plt.plot(X_ls, Y_regr, label=sympy.latex(sympy.Eq(S('y'), L3), mode='inline'))
plt.scatter(X3, Y3, color='red', label='data')
plt.grid()
plt.legend()<jupyter_output><empty_output><jupyter_text>Многочлен Лагранжа с scipy.interpolate.lagrange Пример 5Построим полином Лагранжа для данных Примера 3 с помощью scipy.interpolate.lagrange<jupyter_code>from scipy.interpolate import lagrange
from numpy.polynomial.polynomial import Polynomial
poly = lagrange(X3, Y3)
Polynomial(poly).coef<jupyter_output><empty_output><jupyter_text>Сравним с полученными ранее коэффициентами:<jupyter_code>coeffs = L3.as_coefficients_dict()
res = [round(coeffs[item], 8) for item in [x**k for k in range(len(coeffs))]]
res.reverse()
res<jupyter_output><empty_output><jupyter_text>Квадратичный сплайн Пример 7Аппроксимировать $f(x) = x^3 - 6x^2 + 11 x - 6$ на отрезке $[0, 4]$ квадратичным сплайном с узлами 0, 2, 4.<jupyter_code>def spl2(a0, a1, a2, x):
    return a0 + a1*x + a2*x**2
def dspl2(a0, a1, a2, x, x0):
    return spl2(a0, a1, a2, x).diff(x).subs(x, x0)
def f(x):
    return x**3 - 6*x**2 + 11*x - 6
def df(x, x0):
    return f(x).diff(x).subs(x, x0)
x0 = 0
x1 = 2
x2 = 4
f0, f1, f2 = [f(xi) for xi in (x0, x1, x2)]
df0 = df(x, x0)
a0, a1, a2, x = sympy.symbols('a0:3 x')
sys1 = [sympy.Eq(dspl2(a0, a1, a2, x, x0), df0), sympy.Eq(spl2(a0, a1, a2, x0), f0), sympy.Eq(spl2(a0, a1, a2, x1), f1)]
display(*sys1)
res1 = sympy.solve(sys1)
spline1 = res1[a0] + res1[a1]*x + res1[a2]*x**2
display(res1, spline1)
df1 = spline1.diff(x).subs(x, x1)
sys2 = [sympy.Eq(dspl2(a0, a1, a2, x, x1), df1), 
                   sympy.Eq(spl2(a0, a1, a2, x1), f1), 
                   sympy.Eq(spl2(a0, a1, a2, x2), f2)]
res2 = sympy.solve(sys2)

spline2 = res2[a0] + res2[a1]*x + res2[a2]*x**2
display(*sys2, spline2)
x0, x1, x2 = (0, 2, 4)
X_ls1 = np.linspace(x0, x1)
X_ls2 = np.linspace(x1, x2)
X_ls12 = np.linspace(x0, x2)
Y1 = [spline1.subs(x, item) for item in X_ls1] 
Y2 = [spline2.subs(x, item) for item in X_ls2]          
plt.plot(X_ls1, Y1, color='green')
plt.plot(X_ls2, Y2, color='red')
plt.plot(X_ls12, f(X_ls12), color='black')<jupyter_output><empty_output><jupyter_text>Кубический сплайн$$f_i(x) = a_ix^3 + b_ix^2 + c_ix + d,\quad\left\{\begin{matrix}f_i(x_{i - 1}) = y_{i - 1}\\f_i(x_{i}) = y_{i}\\f_i'(x_{i - 1}) = f'_{i - 1}(x_{i - 1})\\f_i''(x_{i - 1}) = f''_{i - 1}(x_{i - 1})\end{matrix}\right.$$ Пример 8Аппроксимировать $f(x) = \sin(x)$ на отрезке $\left[0, \frac{2\pi}{3}\right]$ квадратичным сплайном с узлами 0, $\frac{\pi}{3}$, $\frac{2\pi}{3}$.<jupyter_code>X = [np.pi*k/3 for k in range(3)]
Y = [np.sin(xk) for xk in X]
plt.scatter(X, Y)<jupyter_output><empty_output><jupyter_text>Вначале определим несколько вспомогательных функций:spl3(a0, a1, a2, a3, x) кубическая парабола с коэффициентами a0, a1, a2, a3 и аргументом xdspl3(a1, a2, a3, x) - первая производная кубического сплайнаd2spl3(a2, a3, x) - вторая производная кубического сплайнаsys_spl(f, x, xi, xi_1, yi, yi_1, d1spl, d2spl) - система уравнений относительно коэффициентов сплайна<jupyter_code>def spl3(a0, a1, a2, a3, x):
    return a0 + a1*x + a2*x**2 + a3*x**3
def dspl3(a1, a2, a3, x):
    return a1 + 2*a2*x + 3*a3*x**2
def d2spl3(a2, a3, x):
    return 2*a2 + 6*a3*x
def sys_spl(f, x, xi, xi_1, yi, yi_1, d1spl, d2spl):
    return [sympy.Eq(spl3(a0, a1, a2, a3, xi_1), yi_1),
            sympy.Eq(spl3(a0, a1, a2, a3, xi), yi),
            sympy.Eq(dspl3(a1, a2, a3, xi_1), d1spl),
            sympy.Eq(d2spl3(a2, a3, xi_1), d2spl)]<jupyter_output><empty_output><jupyter_text>Нам понадобится $\sin(x)$ из numpy для построения графика и из sympy для получения сплайна (аналитически), поэтому опишем функцию fx(x, lib='sympy'), в которой по умолчанию $\sin(x)$ берется из из sympy, но можно изменить значение по умолчанию параметра lib на 'numpy'.<jupyter_code>def fx(x, lib='sympy'):
    if lib == 'sympy':
        return sympy.sin(x)
    if lib == 'numpy':
        return np.sin(x)<jupyter_output><empty_output><jupyter_text>Введем переменные в математическом смысле - символы a0, a1, a2, a3, x.Обозначим xi_1 горизонтальную координату первой точки, xi - второй, yi_1 и yi - соответствующие вертикальные координаты, d1spl и d2spl - значения первой и второй производной на левом конце отрезка, т.е. в xi_1, они для самого левого отрезка равны производным самой функции в этой точке. Для построения сплайна воспользуемся решением с помощью solve соответствующей системы уравнений относительно коэффициентов сплайна.Полученные значения коэффициентов используем для получения выражения (expression) для сплайна на левом отрезке spline3<jupyter_code>a0, a1, a2, a3, x = sympy.symbols('a0:4 x')
xi_1 = X[0]
xi = X[1]
yi_1 = fx(xi_1)
yi = fx(xi)
d1spl = fx(x).diff(x).subs(x, xi_1) 
d2spl = fx(x).diff(x, 2).subs(x, xi_1)
sys3 = sys_spl(fx, x, xi, xi_1, yi, yi_1, d1spl, d2spl)
display(*sys3)
res3 = sympy.solve(sys3)
spline3 = res3[a0] + res3[a1]*x + res3[a2]*x**2 + res3[a3]*x**3
display(spline3)<jupyter_output><empty_output><jupyter_text>Повторим те же действия со следующим отрезком, т.е. $\left[\frac{\pi}{3}, \frac{2\pi}{3}\right]$.Отличие в том, что производные первого и второго порядков считаются не по исходной функции, а по построенному куску сплайна, т.е. считаем производные от  spline3 в точке $\frac{\pi}{3}$.<jupyter_code>xi_1 = X[1]
xi = X[2]
yi_1 = fx(xi_1)
yi = fx(xi)
d1spl = spline3.diff(x).subs(x, xi_1) 
d2spl = spline3.diff(x, 2).subs(x, xi_1)
sys4 = sys_spl(fx, x, xi, xi_1, yi, yi_1, d1spl, d2spl)
display(*sys4)
res4 = sympy.solve(sys4)
spline4 = res4[a0] + res4[a1]*x + res4[a2]*x**2 + res4[a3]*x**3
display(spline4)<jupyter_output><empty_output><jupyter_text>Построим на графике исходную функцию и весь сплайн (из двух кусков)<jupyter_code>x0, x1, x2 = X
X_ls1 = np.linspace(x0, x1)
X_ls2 = np.linspace(x1, x2)
X_ls12 = np.linspace(x0, x2)
Y1 = [spline3.subs(x, item) for item in X_ls1] 
Y2 = [spline4.subs(x, item) for item in X_ls2]          
plt.plot(X_ls1, Y1, color='green', label='spline1')
plt.plot(X_ls2, Y2, color='red', label='spline2')
plt.plot(X_ls12, fx(X_ls12, lib='numpy'), color='black', linestyle=':', label='sin(x)')
plt.legend()<jupyter_output><empty_output><jupyter_text>Сплайны scipy.interpolate: interp1d, splrep, InterpolatedUnivariateSplinehttps://docs.scipy.org/doc/scipy/reference/interpolate.htmlmodule-scipy.interpolate Кубический сплайн:https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.CubicSpline.htmlscipy.interpolate.CubicSpline  interp1dhttps://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.htmlscipy.interpolate.interp1d  Полиномы Бернштейнаhttps://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.BPoly.htmlscipy.interpolate.BPolyТеперь будем строить сплайны средствами scipy.interpolate, вначале подключим необходимые функции:<jupyter_code>from scipy.interpolate import interp1d, splrep, splev, InterpolatedUnivariateSpline, BPoly, CubicSpline<jupyter_output><empty_output><jupyter_text>interp1dПараметры $x$ и  $y$ - координаты узлов сплайна, kind - str или int, необязательный параметр, определяющий тип сплайна, по умолчанию линейный (‘linear’), можно сделать ‘nearest’, ‘nearest-up’, ‘zero’, ‘slinear’, ‘quadratic’, ‘cubic’, ‘previous’ или ‘next’. ‘zero’, ‘slinear’, ‘quadratic’ и ‘cubic’ обозначают интерполяцию нулевого, первого, второго или третьего порядка,‘previous’ и ‘next’ просто возвращают предыдущее или следующее значение функции‘nearest-up’ и ‘nearest’ различаются обработкой полуцелых чисел, таких как 0.5, 1.5. Пример 9Аппроксимировать $f(x) = \sin(x)$ на отрезке $\left[0, 2\pi\right]$ квадратичным сплайном с узлами $\frac{\pi k}{4}$, $k = 0, ..., 8$.<jupyter_code>X = [np.pi*k/4 for k in range(9)]
Y = [np.sin(xk) for xk in X]
plt.scatter(X, Y, color='red')
spl1 = interp1d(X, Y)
spl2 = interp1d(X, Y, kind='cubic')
xs = np.linspace(X[0], X[-1], 1000)
plt.plot(xs, np.sin(xs), 'c-', xs, spl1(xs), 'g:', xs, spl2(xs), 'm--', lw=3)<jupyter_output><empty_output><jupyter_text>Еще один способ построения сплайна: splrep функция, которая возвращает представление  B-spline кривойsplev вычисляет значения B-spline или его производныхЕсть еще полезные функции:sproot находит корни кубического B-splinesplint вычисляет определенный интеграл  B-spline на отрезкеspalde вычисляет все производные B-splinesplrep, splev работают на основе FITPACK (написан на Фортране) Пример 10. Построим сплайн для Примера 7 с помощью splrep и splev.<jupyter_code>X = [np.pi*k/4 for k in range(9)]
Y = [np.sin(xk) for xk in X]
tck = splrep(X, Y, s=0)
xnew = np.linspace(X[0], X[-1], 1000)
ynew = splev(xnew, tck, der=0)
plt.plot(X, Y, 'ro', xnew, ynew, 'c--', xnew, np.sin(xnew), 'k:', lw=3)
plt.legend(['Linear', 'Cubic Spline', 'sin(x)'])<jupyter_output><empty_output><jupyter_text>InterpolatedUnivariateSplineтоже работает на FITPACK. Пример 11. Построим сплайн для Примера 7 с помощью InterpolatedUnivariateSpline (объектно-ориентированная обертка для функций FITPACK).<jupyter_code>X = [np.pi*k/4 for k in range(9)]
Y = [np.sin(xk) for xk in X]
spl = InterpolatedUnivariateSpline(X, Y)
xs = np.linspace(X[0], X[-1], 1000)
plt.scatter(X, Y, color='red')
plt.plot(xs, np.sin(xs), 'c-', xs, spl(xs), 'g--', lw=3)<jupyter_output><empty_output><jupyter_text>Пример 12Построим сплайн для Примера 7 с помощью CubicSpline<jupyter_code>X = [np.pi*k/4 for k in range(9)]
Y = [np.sin(xk) for xk in X]
cs = CubicSpline(X, Y)
xs = np.linspace(X[0], X[-1], 1000)
plt.scatter(X, Y, color='red')
plt.plot(xs, np.sin(xs), 'c-', xs, cs(xs), 'g--', lw=3)<jupyter_output><empty_output><jupyter_text>Кривые Безье Линейная кривая Безье$P_0$ и $P_1$ две различные точки, тогда линейная кривая Безье - отрезок, соединяющий эти точки:$$B(t) = P_0 + t(P_1 - P_0) = (1 - t)P_0 + tP_1,\quad 1\le t \le 1 $$ Квадратичная кривая Безье$P_0$, $P_1$ и $P_2$ три различные точки, тогда квадратичная кривая Безье - отрезок, соединяющий эти точки:$$B(t) = (1 - t)^2P_0 + 2(1 - t)tP_1 + t^2P_2 = P_1 + (1 - t)^2(P_0 - P_1) + t^2(P_2 - P_1),\quad 1\le t \le 1 $$ Кубическая кривая Безье$P_0$, $P_1$, $P_2$ и $P_3$ четыре различные точки, тогда кубическая кривая Безье - отрезок, соединяющий эти точки:$$B(t) = (1 - t)^3P_0 + 3(1 - t)^2tP_1 + 3(1 - t)t^2P_2 +t^3P_3,\quad 1\le t \le 1 $$ Общий вид кривой Безье$P_0$, $P_1$, ...  и $P_n$ - $n$ различных точек, тогда линейная кривая Безье - отрезок, соединяющий эти точки:$$B(t) = \sum_{k=0}^n C_n^k(1 - t)^{n - k}t^kP_k,\quad 1\le t \le 1 $$https://bezier.readthedocs.io/en/stable/python/reference/bezier.curve.html Пример 13Построим кривую Безье по точкам $P_1(-1, 3)$, $P_2(0, 4)$, $P_3(3, 2)$.Для каждой точки построим полином Бернштейна $C_n^k(1 - t)^{n - k}t^k$, $1\le t \le 1$ и составим сумму произведений соответствующих координат точки на полиномы  Бернштейна.<jupyter_code>from scipy.special import comb
def my_B(X, Y):
    n = len(X)
    def my_B_x(t):
        return sum([X[k]*comb(n - 1, k)*(1 - t)**(n - k - 1)*t**k for k in range(n)])
    def my_B_y(t):
        return sum([Y[k]*comb(n - 1, k)*(1 - t)**(n - k - 1)*t**k for k in range(n)])
    return (my_B_x, my_B_y)      
P1 = (-1, 3)
P2 = (0, 4)
P3 = (3, 2)
X = []
Y = []
for point in (P1, P2, P3):
    X.append(point[0])
    Y.append(point[1])
B_x, B_y = my_B(X, Y)
X1 = np.array(X).reshape(3, 1)
Y1 = np.array(Y).reshape(3, 1)
x = [0, 1]
bpX = BPoly(X1, x)
bpY = BPoly(Y1, x)
t_linspace = np.linspace(0, 1)
plt.plot(bpX(t_linspace), bpY(t_linspace), 'c-', B_x(t_linspace), B_y(t_linspace), 'r:', lw=3)
plt.scatter(X, Y)<jupyter_output><empty_output><jupyter_text>Пример 14Построим кривую Безье с помощью модуля bezier:https://bezier.readthedocs.io/en/stable/python/reference/bezier.curve.html<jupyter_code>import bezier
nodes = np.asfortranarray([
...     [0.0, 0.625, 1.0],
...     [0.0, 0.5  , 0.5],
... ])
>>> curve = bezier.Curve(nodes, degree=2)
>>> curve<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17334662678630774, 'fraction_numerical': 0.04885343968095713, 'mean_word_length': 4.603724394785847, 'pattern_counts': {'":': 0, '<': 92, '<?xml version=': 0, '>': 98, 'https://': 6, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '23', 'dir': 'jupyter-structured-clean-dedup', 'id': '78308', 'license': 'MIT', 'n_tokens_mistral': 7425, 'n_tokens_neox': 7672, 'n_words': 1566, 'path': 'Applied_Alg_sem_5_Interpolation.ipynb', 'repo_name': 'GalinaZh/Appl_alg2021', 'size': '15045'}