Announce

PukiWiki contents have been moved into SONOTS Plugin (20070703)

Code: Newton Interpolation

Table of Contents

Abstract

In the mathematical field of numerical analysis, a Newton polynomial, named after its inventor Isaac Newton, is the interpolation polynomial for a given set of data points in the Newton form. The Newton polynomial is sometimes called Newton's divided differences interpolation polynomial because the coefficients of the polynomial are calculated using divided differences.

quoted from wikipedia.

run_newtoninter_chebyshev2.png

First Edition: 2007/02. Last Modified: 2007/02.
Tag: Scientific NumericalAnalysis Interpolation Newton

Source Code

newtoninter.m

function [f, a, d] = newtoninter(x, y, p)
% Newton interpolation
%
%  [f a d] = newtoninter(x, y, p)
%
% Input arguments ([]s are optional):
%  x   (vector) of size 1xN which contains the interpolation nodes.
%  y   (vector) of size 1xN which contains the function values at x
%  p   (vector) of size 1xP which contains points to be interpolated. 
%
% Output arguments ([]s are optional):
%  f   (vector) of size 1xP. The result of interpolation respect to p.
%  [a] (vector) of size 1xN which is leading coefficients genereated by 
%       divided difference method.
%  [d] (matrix) of size NxN (triangular) which is the result of the 
%       divided difference method
%
% Example
% >> x=[1,2,4,7,8]
% >> y=exp(x);
% >> [f a d]= newtoninter(x, y, 5)
%
% Author: Naotoshi Seo <sonots(at)umd.edu>
% Date  : Feb 2007

n = length(x);
d(:,1)=y';
for j=2:n
    for i=j:n
        d(i,j)= ( d(i-1,j-1)-d(i,j-1)) / (x(i-j+1)-x(i));
    end
end
a = diag(d)';

Df(1,:) = repmat(1, size(p));
c(1,:) = repmat(a(1), size(p));
for j = 2 : n
   Df(j,:)=(p - x(j-1)) .* Df(j-1,:);
   c(j,:) = a(j) .* Df(j,:);
end
f=sum(c);

Experiments

Equally Spaced Nodes

Code

function run_newtoninter_equallyspaced
% Date: Feb 2007
% Author: Naotoshi Seo <sonots(at)umd.edu>

% sampling points: eqaully spaced
a = -1; b = 1;
n = 10;
h = (b-a)/(n-1);
x = a + (0:n-1)*h;

% interpolation points: eqaully spaced
n = 50;
h = (b-a)/(n-1);
p = a + (0:n-1)*h;

% (1) f(x) = exp(x)
y = exp(x);
truth = exp(p);
f = newtoninter(x,y,p);

hold on;
plot(x,y,'ob');
plot(p,truth,'-b');
plot(p,f,'-r');
title('Newton Interpolation: Uniform spaced nodes');
legend('samples', 'truth', 'interpolation');

% (2) f(x) = 1 / (1 + 25 * x^2)
figure;
y = 1 ./ (1 + 25 .* x.^2);
truth = 1 ./ (1 + 25 .* p.^2);
f= newtoninter(x,y,p);

hold on;
plot(x,y,'ob');
plot(p,truth,'-b');
plot(p,f,'-r');
title('Newton Interpolation: Uniform spaced nodes');
legend('samples', 'truth', 'interpolation');

Results

run_newtoninter_equallyspaced1.png

run_newtoninter_equallyspaced2.png

Chebyshev nodes

Code

function run_newtoninter_chebyshev
% Date: Feb 2007
% Author: Naotoshi Seo <sonots(at)umd.edu>

% sampling points: chebyshev nodes [-1 1]
a = -1; b = 1;
n = 10;
x = cos( (2 * (1:n) - 1) / (2 * n) * pi);

% interpolation points: eqaully spaced
n = 50;
h = (b-a)/(n-1);
p = a + (0:n-1)*h;

% (1) f(x) = exp(x)
y = exp(x);
truth = exp(p);
f = newtoninter(x,y,p);

hold on;
plot(x,y,'ob');
plot(p,truth,'-b');
plot(p,f,'-r');
title('Newton Interpolation: Chebyshev nodes');
legend('samples', 'truth', 'interpolation');

% (2) f(x) = 1 / (1 + 25 * x^2)
figure;
y = 1 ./ (1 + 25 .* x.^2);
truth = 1 ./ (1 + 25 .* p.^2);
f= newtoninter(x,y,p);

hold on;
plot(x,y,'ob');
plot(p,truth,'-b');
plot(p,f,'-r');
title('Newton Interpolation: Chebyshev nodes');
legend('samples', 'truth', 'interpolation');

Result

run_newtoninter_chebyshev1.png

run_newtoninter_chebyshev2.png

Reference