 ## Announce

Puki contents have been moved into SONOTS Plugin (20070703)

## 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. ## 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  ### 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  