Newton’s Interpolation in MATLAB (Forward and Backward)

Named after Sir Isaac Newton, Newton’s Interpolation is a popular polynomial interpolating technique of numerical analysis and mathematics. Here, the coefficients of polynomials are calculated by using divided difference, so this method of interpolation is also known as Newton’s divided difference interpolation polynomial. Newton polynomial interpolation consists of Newton’s forward difference formula and Newton’s backward difference formula.

In this tutorial, we’re going to write Matlab programs for Newton’s forward interpolation as well as Newton’s backward interpolation, going through the mathematical derivation of the interpolation technique in general. You can also check out our earlier tutorials where we discussed C programs for Newton’s forward interpolation and Newton’s divided difference.

Derivation of Newton’s Interpolation:

Let us consider a set of k+1 data points (x₀, y₀) , (x₁, y₁), (x₂, y₂), . … . .(x_k, y_k).

where, every x_j is unique.

The given interpolation of polynomial in Newton’s form can be expressed in linear combination of Newton basis polynomial as follows:

Again, the Newton basis Polynomials can be defined as:

which is valid for j > 0 and n₀(x) ≡ 1

Similarly, the coefficients are defined as:

a_j := [y₀, y₁, y₂, .. . . . . .y_j]

where, [y₀, y₁, y₂, .. . . . . .y_j] is the notation of divided difference.

Now, Newton’s interpolation or polynomial can be expressed as:

N(x) = [y₀] + [y₀, y₁] (x –x₀) + . . . + [y₀, . . .y_k]( x – x₀) ( x – x₁) . . . ( x – x_k-1)

This form of Newton’s polynomial can be simplified by arranging x_o, x₁, x₂, … x_k in consecutively equal space.

For simplicity, use the notation h = x_i+1 – x_i

where, i = 0, 1, . . . .k- 1 and x = x₀ + sh

Again, the difference x – xi is written as ( s – i )h.

Finally, the above expression for Newton polynomial above becomes:

N(x) = [y₀] + [ y₀, y₁]sh + . . . + [y₀, y₁, . . . .y_k] s(s – 1) . . . (s – k +1) h^k

which is called Newton’s forward difference formula.

If the nodes are recorded from backward as x_k, x_k-1, . . . , x₀

The Newton’s Polynomial can be expresses as:

N(x) = [y_k] + [y_k, y_k-1] (x – x_k) + . . . + [y_k, y_k-1, . . . y₀] (x – x_k) ( x – x_k-1) . . . (x – x₁)

If the xk, x_k-1, . . .  . x0 are equally spaced and x = x_k + sh and x_i = x_k– ( k –i)h for i = 0, 1, . . k

Now, Newton’s polynomial becomes:

N(x) = [y_k] + [ y_k, y_k-1]sh + . . . + [ y_k, y_k-1, . . . .y₀]s(s+1) . . . (s+ k – 1) h^k

which is called Newton’s Backward Difference Formula.

Newton’s Interpolation in MATLAB:

Here are two different MATLAB codes for Newton’s forward as well as backward interpolation, written on the basis of aforementioned derivation cum formula.

Forward Interpolation using MATLAB:

%Newton’s Forward Difference Formula MATLAB Program
x=[0 2 4 7 10 12]; % inputting values of x 
fx=[20 20 12 7 6 6]; % inputting values of y
dt=zeros(6,10); % function 
for i=1:6 dt(i,1)=x(i);% for loop 
dt(i,2)=fx(i); % calling function 
end 
n=5; % number of iterations 
for j=3:10
for i=1:n
dt(i,j)=dt(i+1,j-1)-dt(i,j-1)
end
n=n-1;
end 
h=x(2)-x(1) % finding the value of h
xp=1.5; % defining the value of xp 
for i=1:5
q=(xp-x(i))/h; % calculating number of intervals 
if (q>0&&q<1)
p=q;
 end
 end
 p
l=xp-(p*h)
 for i=1:5 
 if(l==x(i))
 r=i;
 end
 end % calculating different value of y
 f0=fx(r);
 f01=dt(r,3); 
 f02=dt(r,(3+1));
 f03=dt((r),(3+2));
 f04=dt((r),(3+3));
% using the forward interpolation formula 
 fp=(f0)+((p*f01)+(p*(p-1)*f02)/(2)) + ((p*(p-1)*(p-2)*f03)/(6))+((p*(p-1)*(p-2)*(p-3)*f04)/(24))

The above source code for Newton’s Interpolation using forward difference formula doesn’t need any input value. All the input values required for the interpolation are embedded within the source code. The values of x and y used in above source code are 0 2 4 7 10 12 and 20 20 12 7 6 6 respectively. The sample output of the program is given below:

Backward Interpolation using MATLAB:

%Newton’s Backward Difference Formula MATLAB Program
x=[0 8 16 24 32 40]; % inputting the values of x 
fx=[14.621 11.843 9.870 8.418 7.305 6.413]; % inputting the value of y 
dt=zeros(6,7); % declaring function 
for i=1:6 % stating loop 
dt(i,1)=x(i); 
dt(i,2)=fx(i);
end
n=5; 
for j=3:7
 for i=1:n % using for loop 
 dt(i,j)=dt(i+1,j-1)-dt(i,j-1) % defining dt 
 end n=n-1;
 end
h=x(2)-x(1) % finding the value of h
xp=27; % defining xp
for i=1:6
q=(xp-x(i))/h;
if (q>0&&q<1)
 p=q;
 end
 end
 p
 l=xp-(p*h)
 for i=1:6
 if(l==x(i))
 r=i;
 end
 end
% finding different value of y
 f0=fx(r);
 f01=dt((r-1),3);
 f02=dt((r-2),(3+1));
f03=dt((r-3),(3+2));
 f04=dt((r-4),(3+3));
%using backward difference formula 
 fp=(f0)+((p*f01)+(p*(p+1)*f02)/(2)) + ((p*(p+1)*(p+2)*f03)/(6))

The above source code in MATLAB for Newton’s interpolation by using backward difference formula is to solve the following data:

X: 0 8 16 24 32 40

Y :14.621 11.843 9.870 8.418 7.305 6.413

Since, all the input data are already defined in the source code, the user doesn’t need to input any value to the program while running. The sample output of the above Matlab code is given below:

If you have any questions regarding the above two Newton’s interpolation techniques, their derivation or MATLAB programs, bring them up from the comments section. You can find more Numerical Methods tutorial using MATLAB here.