%
% Newton 插值
%

clear; clc;

X = [0.5, 0.6, 0.4];  % 插值节点
Y = [ -0.6931, -0.5108, -0.9163]; % 函数值
[p,q] = divided_difference(X,Y);   % 计算差商, q 为差商表对角线

x = 0.54;  % 需要插值的点
% 线性插值：取插值节点 0.5, 0.6
y1 = q(1) + q(2)*(x-0.5);
fprintf('线性插值的结果: %.4f\n',y1);

% 二次插值：取插值节点 0.5, 0.6, 0.4
y2 = q(1) + q(2)*(x-0.5) + q(3)*(x-0.5)*(x-0.6);
fprintf('二次插值的结果: %.4f\n',y2);


%
%  计算差商
%
function [p, q] = divided_difference(x,y)
%
% 参数说明：
% 输入参数： x 为节点，y 为函数值（注意：向量 x 与向量 y 的长度必须一致）
% 输出参数：
%    p 包含所有差商，即差商表中的所有值
%    q 为差商表中对角线的值

m = length(x);
x = x(:);
p = zeros(m, m+1);
p(:,1) = x; 
p(:,2) = y(:);
for k = 3 : m+1
    p(k-1:m, k) = diff(p(k-2:m, k-1)) ./ ( x(k-1:m) - x(1:m+2-k) );
end
q = diag(p(1:m,2:m+1));
end