Интерполяционный многочлен Ньютона и Лагранджа

function [L_n P] = lagrange1(x,y,xx)
% x - массив координат узлов
% y - массив значений интерполируемой функции
% xx - массив значений точек интерполяции
L_n = zeros(size(xx));
for k=1:length(x)
    P = ones(length(x),length(xx));
    for i=1:length(x)
        if k~=i
            P(k,:) = P(k,:).*(xx-x(i))./(x(k)-x(i));
        end
    end
    L_n = L_n + y(k)*P(k,:);
end
end

clear, clc
 
x = [0.165 0.171, 0.177];
y = [2.00551 2.82558 2.65583];
xx = linspace(0.16,0.18,100); % массив точек х
x0 = 0.169; % некоторая точка
yy = lagrange1(x,y, xx);
[y0 P0] = lagrange1(x,y, x0) % функция в некоторой точке и полином
plot(x,y,'or',xx,yy,':r',x0,y0,'*b'), grid on
legend('Данные','Интерполяция','Точка', 'location', 'northwest')

function [L_n P0] = lagrange1(x,y,xx)
% x - массив координат узлов
% y - массив значений интерполируемой функции
% xx - массив значений точек интерполяции
L_n = zeros(size(xx));
for k=1:length(x)
    P = ones(size(xx));
    for i=1:length(x)
        if k~=i
            P = P.*(xx-x(i))./(x(k)-x(i));
        end
    end
    P0(k,:) = P;
    L_n = L_n + y(k)*P;
end
end

clear, clc
 
x = [0.165 0.171, 0.177];
y = [2.00551 2.82558 2.65583];
xx = linspace(0.16,0.18,100);
x0 = 0.169;
yy = lagrange1(x,y, xx);
[y0 P0] = lagrange1(x,y, x0)
plot(x,y,'or',xx,yy,':r',x0,y0,'*b'), grid on
legend('Данные','Интерполяция','Точка', 'location', 'northwest')
 
y*P0 % проверка

clear, clc
 
x = [0.165 0.171, 0.177];
y = [2.00551 2.82558 2.65583];
x0 = 0.169;
y0 = lagrange1(x,y,x0);
syms X
P1 = ((X-x(2))*(X-x(3)))/((x(1)-x(2))*(x(1)-x(3)));
P2 = ((X-x(1))*(X-x(3)))/((x(2)-x(1))*(x(2)-x(3)));
P3 = ((X-x(1))*(X-x(2)))/((x(3)-x(1))*(x(3)-x(2)));
Y = y(1)*P1 + y(2)*P2 + y(3)*P3;
 
pretty(collect(Y,X))
Yx = sym2poly(Y)
 
y0 
polyval(Yx,x0)

clear, clc
% координаты узлов
x = [0.165 0.171, 0.177];
y = [2.00551 2.82558 2.65583];
% некая точка
x0 = 0.169;
y0 = lagrange1(x,y,x0);
syms X
for i = 1:length(x)
    u = x;
    u(i)=[]; % выкалываем i-й элемент
    P(i,1) = prod(X-u)/prod(x(i)-u); % рассчитываем i-й полином Лагранжа
end
Y = y*P % сумма произведений y на 
 
pretty(collect(Y,X)) % собираем Х
Yx = sym2poly(Y) % вытягиваем коэфф-ты многочлена
 
y0 % расчитанная ранее функция, для стравнения
polyval(Yx,x0) % подставляем значение х0 в полученный полином