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"
using Plots
f = x-> 1./(1+x.^2) # create a function f
plot(f,-5,5,legend=false)
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
"""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.])
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)
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