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