Are all sets of interpolation points equivalent?

In this demo, we'll look at two ways of interpolating a function $f(x)$. The function we use is a simple one: $$ f(x) = 1/(1+x^2) $$ This is also called "Runge's function"

In [2]:
using Plots
f = x-> 1./(1+x.^2) # create a function f
plot(f,-5,5,legend=false)
[Plots.jl] Initializing backend: pyplot
WARNING: No working GUI backend found for matplotlib.
Out[2]:
In [3]:
N = 7
xx1 = linspace(-5,5,N)
xx2 = 5*cos(pi*(0:(N-1))./(N-1))
plot(f,-5,5,border=false,legend=false,ticks=false)
scatter!(xx1,-0.05*ones(N)) # point set 1
scatter!(xx2,-0.1*ones(N)) # point set 2
Out[3]:
In [4]:
"""fit the coefficients of a polynomial interpolant
c = polyfit(x,y,n) fits a degree n polynomial to the
data poly(x,c)"""
function polyfit(x,y,n)
    m = length(x) # datapoints
    Z = zeros(m,n+1)
    for i=1:m # Should fill this out columnwise.
        xi = 1.
        for j=1:n+1
            Z[i,j] = xi
            xi *= x[i]
        end
    end
    return Z\y
end
""" Evaluate a polynomial vs. Horner's rule """
function polyval(c,x)
    # based on Matlab's Horner's rule 
    m = length(x)
    y = zeros(m)
    
    if length(c) > 0
        y[:] = c[end]
    end
    for i=2:length(c)
        y = x .* y + c[end-i+1]
    end
    return y
end
@show polyfit([1,2,3.],[1,4.,9],2)
@show polyval([0.,0.,1.],[0.5,4.])
polyfit([1,2,3.0],[1,4.0,9],2) = [0.0,0.0,1.0]
polyval([0.0,0.0,1.0],[0.5,4.0]) = [0.25,16.0]
Out[4]:
2-element Array{Float64,1}:
  0.25
 16.0 
In [5]:
N = 15
xx1 = linspace(-5,5,N)
xx2 = 5*cos(pi*(0:(N-1))./(N-1))

# build two interpolants
c1 = polyfit(xx1,f(xx1),N-1)
c2 = polyfit(xx2,f(xx2),N-1)

xxall = linspace(-5,5,1000)
Y = [f(xxall) polyval(c1,xxall) polyval(c2,xxall)]
l = @layout([a;b;c])
plot(xxall,Y,layout=l)
Out[5]:
In [6]:
using Interact
@manipulate for N=3:21
    xx1 = linspace(-5,5,N)
    xx2 = 5*cos(pi*(0:(N-1))./(N-1))

    # build two interpolants
    c1 = polyfit(xx1,f(xx1),N-1)
    c2 = polyfit(xx2,f(xx2),N-1)

    xxall = linspace(-5,5,1000)
    Y = [f(xxall) polyval(c1,xxall) polyval(c2,xxall)]
    plot(xxall,Y)
    scatter!(xx1,f(xx1))
    scatter!(xx2,f(xx2))
end
Out[6]:
In [ ]: