In [1]:
using Plots
In [8]:
""" 
Evaluate the polynomial interpolant in Lagrange form.

xi is a vector of length n+1
fi are the values f(xi) for each point xi

x is the value where we compute
p(x)
which is the degree n+1 polynomial interpolating f at each point xi

The values xi must be distinct

Note that x can be a matrix or vector
"""
function lagrange_interp(x, xi, fi)
    f = zeros(size(x))
    for i=1:length(fi)
        elli = ones(size(x))
        for j=1:length(xi)
            if i==j continue; end
            elli = elli.*(x-xi[j])/(xi[i] - xi[j])
        end
        f = f + fi[i]*elli
    end
    return f
end

# The polynomial
fi = [2,3,6];
xi = [1,2,3];
xx = linspace(0,4,100);

plot(xx,lagrange_interp(xx,xi,fi))
Out[8]:
In [5]:
function compute_barycentric_weights(xi)
    ws = zeros(size(xi))
    for i=1:length(xi)
        w = 1
        for j=1:length(xi)
            if j==i continue; end
            w = w*1./(xi[i] - xi[j]);
        end
        ws[i] = w;
    end
    return ws
end  

function barycentric_interp_with_weights(x,xi,fi,wi)
    fnumer = zeros(size(x))
    fdenom = zeros(size(x))
    for i=1:length(xi)
        fnumer += fi[i].*wi[i]./(x - xi[i])
        fdenom += wi[i]./(x - xi[i])
    end
    f = fnumer./fdenom
end

function barycentric_interp(x,xi,fi)
    return barycentric_interp_with_weights(x,xi,fi,compute_barycentric_weights(xi))
end

# The polynomial
fi = [2,3,6];
xi = [1,2,3];
xx = linspace(0,4,100);

plot(xx,barycentric_interp(xx,xi,fi))
Out[5]:
In [6]:
# show Trefethen's demo
fun(x) = abs(x) + .5.*x- x.^2;
n = 10
xi = cos(pi*(0:n)/n);
fi = fun(xi)
xx = linspace(-1,1,5000)
plot(xx,barycentric_interp(xx,xi,fi))
scatter!(xi,fi)
Out[6]:
In [7]:
using Interact
using Plots
gadfly()
INFO: Recompiling stale cache file /Users/dgleich/.julia/lib/v0.4/Reactive.ji for module Reactive.
Out[7]:
Plots.GadflyPackage()
In [8]:
fun(x) = abs(x) + .5.*x- x.^2;

@manipulate for n=10:10:1000
    xi = cos(pi*(0:n)/n);
    fi = fun(xi)
    xx = linspace(-1,1,5000)
    plot(xx,barycentric_interp(xx,xi,fi))
    scatter!(xi,fi,markersize=1.).o
end
[Plots.jl] Initializing backend: gadfly
Out[8]:
-3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 -3.0 -2.9 -2.8 -2.7 -2.6 -2.5 -2.4 -2.3 -2.2 -2.1 -2.0 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 -4 -2 0 2 4 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 y1 y2 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 -2.00 -1.95 -1.90 -1.85 -1.80 -1.75 -1.70 -1.65 -1.60 -1.55 -1.50 -1.45 -1.40 -1.35 -1.30 -1.25 -1.20 -1.15 -1.10 -1.05 -1.00 -0.95 -0.90 -0.85 -0.80 -0.75 -0.70 -0.65 -0.60 -0.55 -0.50 -0.45 -0.40 -0.35 -0.30 -0.25 -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 1.55 1.60 1.65 1.70 1.75 1.80 1.85 1.90 1.95 2.00 2.05 2.10 2.15 2.20 2.25 2.30 2.35 2.40 2.45 2.50 -2 0 2 4 -2.0 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5

In [9]:
fun(x) = abs(x) + .5.*x- x.^2;

@manipulate for n=5:5:50
    xi = collect(-1.:(2/n):1.)
    fi = fun(xi)
    xx = linspace(-1,1,5000)
    plot(xx,barycentric_interp(xx,xi,fi))
    scatter!(xi,fi,markersize=1.).o
    ylims!(-1.,1.).o
end
Out[9]:
-3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 -3.0 -2.9 -2.8 -2.7 -2.6 -2.5 -2.4 -2.3 -2.2 -2.1 -2.0 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 -4 -2 0 2 4 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 y1 y2 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 -3.0 -2.9 -2.8 -2.7 -2.6 -2.5 -2.4 -2.3 -2.2 -2.1 -2.0 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 -4 -2 0 2 4 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0
In [ ]: