Newton Forward difference method

% Confirmed case Sample data points (x, y) 
data = [ 1, 88; 2, 49; 3, 47; 4, 8; 5, 34; 6, 762; 7, 98; 8, 40];
% Extract x and y values 
x = data(:, 1); 
y = data(:, 2); 
n = length(x);
% Calculate forward differences
forward_diff = zeros(n, n); 
forward_diff(:, 1) = y;
for j = 2:n 
    for i = 1:n-j+1 
        forward_diff(i, j) = forward_diff(i+1, j-1) - forward_diff(i, j-1); 
    end 
end
% Construct the Newton Forward Difference polynomial 
syms X; 
P = y(1); 
for j = 1:n-1 
    term = forward_diff(1, j); 
    for i = 1:j 
        term = term * (X - x(i)); 
    end 
    P = P + term / factorial(j); 
end
disp('Newton Forward Difference Polynomial:'); 

Newton Forward Difference Polynomial:

disp(simplify(P));

Result Newton Forward Difference Polynomial: - (725*X^7)/1008 + (1651*X^6)/80 - (173233*X^5)/720 + (70651*X^4)/48 - (91321*X^3)/18 + (146407*X^2)/15 - (1996741*X)/210 + 3658

The problem now is that when I substitute each value of x, I don't get the corresponding value of y as stated in the table. I need help on how to get the correct polynomial equation.

Answer 1

[Torsten]

Three coding errors:

% Confirmed case Sample data points (x, y) 
data = [ 1, 88; 2, 49; 3, 47; 4, 8; 5, 34; 6, 762; 7, 98; 8, 40];
% Extract x and y values 
x = data(:, 1); 
y = data(:, 2); 
n = length(x);
% Calculate forward differences
forward_diff = zeros(n, n); 
forward_diff(:, 1) = y;
for j = 2:n 
    for i = 1:n-j+1 
        %forward_diff(i, j) = (forward_diff(i+1, j-1) - forward_diff(i, j-1))
        forward_diff(i, j) = (forward_diff(i+1, j-1) - forward_diff(i, j-1)) / (x(i+j-1) - x(i)); 
    end 
end
% Construct the Newton Forward Difference polynomial 
syms X; 
P = y(1); 
for j = 1:n-1 
    %term = forward_diff(1, j); 
    term = forward_diff(1, j+1);
    for i = 1:j 
        term = term * (X - x(i)); 
    end 
    %P = P + term / factorial(j); 
    P = P + term;
end
disp('Newton Forward Difference Polynomial:'); 

Newton Forward Difference Polynomial:

disp(simplify(P));
double(subs(P,X,x))

ans = 8×1

88 49 47 8 34 762 98 40

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

For all data series with deltax = 1 (as the one above), there is one error in your code. But note that this is a special case and the above code holds for arbitrary deltax.

% Confirmed case Sample data points (x, y) 
data = [ 1, 88; 2, 49; 3, 47; 4, 8; 5, 34; 6, 762; 7, 98; 8, 40];
% Extract x and y values 
x = data(:, 1); 
y = data(:, 2); 
n = length(x);
% Calculate forward differences
forward_diff = zeros(n, n); 
forward_diff(:, 1) = y;
for j = 2:n 
    for i = 1:n-j+1 
        forward_diff(i, j) = forward_diff(i+1, j-1) - forward_diff(i, j-1);
    end 
end
% Construct the Newton Forward Difference polynomial 
syms X; 
P = y(1); 
for j = 1:n-1 
    %term = forward_diff(1, j); 
    term = forward_diff(1, j+1);
    for i = 1:j 
        term = term * (X - x(i)); 
    end 
    P = P + term / factorial(j); 
end
disp('Newton Forward Difference Polynomial:'); 

Newton Forward Difference Polynomial:

disp(simplify(P));
double(subs(P,X,x))

ans = 8×1

88 49 47 8 34 762 98 40

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>