Untitled

function dd = divided_differences(x, y)
    n = length(x);
    dd = zeros(n, n);
    dd(:,1) = y';  // Pierwsza kolumna to wartości y

    for j = 2:n
        for i = 1:n-j+1
            dd(i,j) = (dd(i+1,j-1) - dd(i,j-1)) / (x(i+j-1) - x(i));
        end
    end
endfunction

// Funkcja do interpolacji Newtona
function yi = newton_interpolation(x, dd, xi)
    n = length(x);
    yi = zeros(size(xi));

    for k = 1:length(xi)
        yi(k) = dd(1,1);  
        p = 1;
        for j = 1:n-1
            p = p * (xi(k) - x(j));
            yi(k) = yi(k) + dd(1,j+1) * p;
        end
    end
endfunction

// Funkcja do obliczania pochodnej funkcji interpolującej w punkcie xi
function dyi = newton_interpolation_derivative(x, dd, xi)
    n = length(x);
    dyi = zeros(size(xi));
    for k = 1:length(xi)
        p = 1;
        dp = 0;
        for j = 1:n-1
            dp = dp * (xi(k) - x(j)) + p;
            p = p * (xi(k) - x(j));
            dyi(k) = dyi(k) + dd(1,j+1) * dp;
        end
    end
endfunction

// Dane
x = [-4 -1 0 -1 5 6 8 10];
y = [3 1 8 -1 -5 -4 3 8];

dd = divided_differences(x, y);

xi = linspace(min(x), max(x), 100); 
yi = newton_interpolation(x, dd, xi);

clf; 
plot(x, y, 'ro', 'linewidth', 2);  
plot(xi, yi, '-b', 'linewidth', 2);

xtitle('Interpolacja Newtona', 'X', 'Y');

// Obliczanie pochodnej w punkcie x = 2
x_value = 0;
dyi = newton_interpolation_derivative(x, dd, [x_value]);

disp("Pochodna w punkcie x = 0:");
disp(dyi);