trying to create a function for simple Newton polynomial interpolation

I am trying to create a function for simple Newton polynomial interpolation.

This is the recurrance relation I am trying to create on matlab, where the data points (x_i,y_i) are specified by the input vectors x and y, and the points at which I want to evaluate the polynomial at are given by the input vector t. I then want to store the evaluated points in a vector p and display this vector.

I have been inputting two vectors of the same size x and y and then letting the vector t = x so that if everything went right I should get the vector y back out, but I don't and I'm not sure why? Any Help.

nt = length ( t );

%Number of data points
n = length ( x );


%Initial polynomial approximation
p = y(1);


for k = 2:n
 %initialise product in the recursive formula
 
for j = 1:k-1 
p = p + (y(k)-p(x(k)))*(t-x(j))/(x(k)-x(j));
end

end

p = p (1:nt);

for l = 1:nt
p(nt) = p(t(nt));

end

display(p);

end

Answer 1

[Shashi Kiran]

Hi @Ole,

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

Evaluated polynomial values at points in t:
     1    49     9    16

Hope this helps.​