In [1]:
using Plots

In [3]:
"""
Compute the divided difference tables to get the Newton coefficients.
This function stores the whole table for pedagogical reasons, not
a good general purpose implementation.

This function will return the coeffients on the diagonal.
"""
function compute_newton_coeffs(xi,fi)
npts = length(xi)
T = zeros(npts,npts)
T[:,1] = fi # the first column becomes fi.
for j=2:npts
for i=j:npts
# the value in the table is
# made by looking left (i,j-1) and left-up (i-1,j-1)
# divided by the region x_i - ...
T[i,j] = (T[i,j-1] - T[i-1,j-1])/(xi[i] - xi[i-j+1])
end
end
return T
end

# This is the example from Gautschi, page 96
# f(x) = 3+x^2
xi = [0., 1., 2., 4.]
fi = [3., 4., 7., 19.]
compute_newton_coeffs(xi,fi)

Out[3]:
4x4 Array{Float64,2}:
3.0  0.0  0.0  0.0
4.0  1.0  0.0  0.0
7.0  3.0  1.0  0.0
19.0  6.0  1.0  0.0
In [5]:
function newton_interp(xx,xi,fi)
T = compute_newton_coeffs(xi,fi)
xv = ones(size(xx))
f = T[1,1].*xv
for i=2:size(T,1)
xv = xv.*(xx - xi[i-1])
f += T[i,i]*xv
end
f
end

# Reuse the same example
xx = collect(linspace(-1.,5.,100));
plot(xx,newton_interp(xx,xi,fi))

Out[5]:
In [14]:
fun(x) = 3+0.5*x.^2+exp(-x)
xi = [1.;2+collect(linspace(0.,0.00001,3)); 3.; 4.; ]
fi = fun(xi)
@show xi
T = compute_newton_coeffs(xi,fi)

xi = [1.0,2.0,2.000005,2.00001,3.0,4.0]
Out[14]:
6x6 Array{Float64,2}:
3.86788  0.0      0.0        0.0        0.0          0.0
5.13534  1.26746  0.0        0.0        0.0          0.0
5.13534  1.86467  0.597209   0.0        0.0          0.0
5.13535  1.86467  0.567654  -0.0295546  0.0          0.0
7.54979  2.41446  0.549787  -0.017867   0.00584382   0.0
11.0183   3.46853  0.527038  -0.0113743  0.00324635  -0.000865823

In [2]:
LoadError: UndefVarError: compute_newton_coeffs not defined
while loading In[2], in expression starting on line 3