To create a function that performs Newton polynomial interpolation based on the recurrence relation given these steps can be followed.
- Start with the base polynomial
where 
is the first y-value. - This base value will serve as the initial polynomial value at each point in the vector t where we want to evaluate the polynomial.
function p = newton_interpolation_algorithm(x, y, t)
n = length(x);
nt = length(t);
p = zeros(1, nt); % Initialize output vector for evaluated points
for j = 1:nt
p(j) = y(1); % Initialize p(t(j)) as y_0 for each evaluation point t(j)
end
- For each evaluation point
, calculate a product term that captures the influence of all previous 
values: 
- Use this product term to update
with a new term derived from 
- After computing the product term, add the new term to
: 
% Loop for each additional data point
for k = 2:n
for j = 1:nt
% Calculate product term for each t(j)
product_term = 1;
for i = 1:k-1
product_term = product_term * (t(j) - x(i)) / (x(k) - x(i));
end
% Update polynomial value at t(j)
p(j) = p(j) + (y(k) - p(j)) * product_term;
end
end
- Once the recursion is complete, p contains the evaluated polynomial values at each point in t.
disp('Evaluated polynomial values at points in t:');
disp(p);
end
x = [1, 7, 3, 4];
y = [1, 49, 9, 16];
t = x; % Points to evaluate the polynomial (for testing)
p = newton_interpolation_algorithm(x, y, t); % Should return values close to y if correct
Hope this helps.
