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

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]:

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]:
In [ ]: