26) Chebyshev Polynomials and Piecewise Interpolation#
Last time:#
Interpolation using polynomials and function approximation
Lagrange Interpolating Polynomials
Newton polynomials
Choice of points and choice of bases
Runge phenomenon as a manifestation of ill-conditioning
Today:#
Chebyshev polynomials
1.1 Chebyshev polynomials are well-conditionedIntroduction to Piecewise Interpolation
using LinearAlgebra
using Plots
using LaTeXStrings
default(linewidth=4, legendfontsize=12)
function vander(x, k=nothing)
if isnothing(k)
k = length(x)
end
m = length(x)
V = ones(m, k)
for j in 2:k
V[:, j] = V[:, j-1] .* x
end
V
end
function vander_newton(x, abscissa=nothing)
if isnothing(abscissa)
abscissa = x
end
n = length(abscissa)
A = ones(length(x), n)
for i in 2:n
A[:,i] = A[:,i-1] .* (x .- abscissa[i-1])
end
A
end
CosRange(a, b, n) = (a + b)/2 .+ (b - a)/2 * cos.(LinRange(-pi, 0, n))
CosRange (generic function with 1 method)
So far, we have seen different strategies for Interpolation.
Note
Interpolation \(\neq\) Function Approximation.
But can we use interpolation to accurately approximate continuous functions?
Suppose the interpolation data come from samples of a continuous function \(f\) on \([a, b] \in \R\).
Then we’d like the interpolant to be “close to” \(f\) on \([a, b]\).
1. Chebyshev polynomials#
We saw that interpolating a set of data points using the monomial basis leads to the Vandermonde matrix. We noticed that when using the monomial basis for \(P_n\), that is, \({1, x, x^2, \dots, x^n}\), these monomial basis functions become increasingly similar. This means that the Vandermonde matrix columns become nearly linearly-dependent, i.e., it is an ill-conditioned matrix (hence, it is very sensitive to perturbations in the data).
A well-conditioned polynomial basis:#
Let’s construct the Vandermonde matrix, but instead of using the canonical basis for polynomials, we’ll use the Legendre polynomials.
function vander_legendre(x, k=nothing)
if isnothing(k)
k = length(x) # Square by default
end
m = length(x)
Q = ones(m, k)
Q[:, 2] = x
for n in 1:k-2
Q[:, n+2] = ((2*n + 1) * x .* Q[:, n+1] - n * Q[:, n]) / (n + 1) # definition of Legendre's polynomials via the generating function
end
Q
end
vander_legendre (generic function with 2 methods)
vcond(mat, points, nmax) = [cond(mat(points(-1, 1, n))) for n in 2:nmax]
plot([vcond(vander_legendre, CosRange, 20)], yscale=:log10, label = "cond of vander_legendre with cos nodes")
A different set of points: Chebyshev nodes#
CosRange(a, b, n) = (a + b)/2 .+ (b - a)/2 * cos.(LinRange(-pi, 0, n))
plot([vcond(vander, LinRange, 20)], yscale=:log10, size=(1000, 600), label = "vander LinRange")
plot!([vcond(vander, CosRange, 20)], yscale=:log10, label = "vander CosRange")
plot!([vcond(vander_legendre, LinRange, 20)], yscale=:log10, label = "vander_legendre LinRange")
plot!([vcond(vander_legendre, CosRange, 20)], yscale=:log10, label = "vander_legendre CosRange", legend = :topleft)
What’s wrong with ill-conditioning?#
runge(x) = 1 / (1 + 10*x^2)
x_lin = LinRange(-1, 1, 20)
x_cos = CosRange(-1, 1, 20)
y_lin = runge.(x_lin)
y_cos = runge.(x_cos)
plot(runge, xlims=(-1, 1))
scatter!(x_lin, y_lin)
plot(runge, xlims=(-1, 1))
scatter!(x_cos, y_cos)
ourvander = vander_legendre
p = ourvander(x_cos) \ y_cos
scatter(x_cos, y_cos, label = L"(x_i,y_i)")
s = LinRange(-1, 1, 100)
plot!(s, runge.(s), label = "runge")
plot!(s, ourvander(s, length(x_cos)) * p, label = "vander_legendre")
# Let's check conditioning with the Chebyshev nodes:
svdvals(ourvander(s, length(x_cos)) / ourvander(x_cos))
20-element Vector{Float64}:
2.8563807370813317
2.841466860881556
2.8202363866684843
2.782971972925489
2.7474280450222
2.6846513464302872
2.636781366412176
2.545083971162913
2.4861340570680777
2.3619485533713127
2.2916184164905107
2.1320708818800687
2.0460305389351454
1.854734771453593
1.7347507488398968
1.5501994693431675
1.3481998183116586
1.2519073527674072
0.9462546371465372
0.9419175178348708
cond(ourvander(x_cos))
8.502660254099553
Definition of Chebyshev polynomials#
Define
Recall
Yet another basis:#
Let’s construct again the Vandermonde matrix, but instead of using the canonical basis for polynomials, we’ll use the Chebyshev polynomials.
function vander_chebyshev(x, n=nothing)
if isnothing(n)
n = length(x) # Square by default
end
m = length(x)
T = ones(m, n)
if n > 1
T[:, 2] = x
end
for k in 3:n
T[:, k] = 2 * x .* T[:,k-1] - T[:, k-2]
end
T
end
vander_chebyshev (generic function with 2 methods)
x = LinRange(-1, 1, 50)
plot(x, vander_chebyshev(x, 5), title="Chebyshev polynomials")
1.1 Chebyshev polynomials are well-conditioned#
plot_vcond(mat, points) = plot!([
cond(mat(points(-1, 1, n)))
for n in 2:30], label="$mat/$points", marker=:auto, yscale=:log10)
plot(title="Vandermonde condition numbers", legend=:topleft, size=(1000, 600))
plot_vcond(vander, LinRange)
plot_vcond(vander, CosRange)
plot_vcond(vander_chebyshev, LinRange)
plot_vcond(vander_chebyshev, CosRange)
Lagrange interpolating polynomials revisited#
Let’s re-examine the Lagrange polynomials as we vary the points.
When using Chebyshev nodes in a Vandermonde matrix, the corresponding Lagrange basis functions are constructed using these nodes, allowing for a well-conditioned interpolation process.
x_cos = CosRange(-1, 1, 11)
x_lin = LinRange(-1, 1, 100)
A = vander_chebyshev(x_lin, length(x_cos)) /
vander_chebyshev(x_cos)
plot(x_lin, A[:,1:5], size=(800, 600))
scatter!(x_cos[1:5], ones(5), color=:black, label=nothing)
scatter!(x_cos, zero(x_cos), color=:black, label=nothing)
y = abs.(x_cos) # not a smooth function
plot(x_lin, [abs.(x_lin) A * y])
Are there artifacts?
2. Introduction to Piecewise Interpolation#
Main idea: Break the domain into subdomains, apply polynomial interpolation on each subdomain (interpolation points/nodes now called “knots”).
The simplest piecewise interpolation is the piecewise constant interpolation.
function interp_nearest(x, s)
A = zeros(length(s), length(x))
for (i, t) in enumerate(s)
loc = nothing
dist = Inf
for (j, u) in enumerate(x)
if abs(t - u) < dist
loc = j
dist = abs(t - u)
end
end
A[i, loc] = 1
end
A
end
interp_nearest(LinRange(-1, 1, 3), LinRange(0, 1, 4))
4×3 Matrix{Float64}:
0.0 1.0 0.0
0.0 1.0 0.0
0.0 0.0 1.0
0.0 0.0 1.0
A = interp_nearest(x_cos, x_lin)
plot(x_lin, A[:, 1:5])
scatter!(x_cos, ones(length(x_cos)))
Are there any visual artifacts?
x = LinRange(-1, 1, 40)
A = interp_nearest(x, x_lin)
cond(A)
1.224744871391589
plot(x_lin, A * runge.(x))
plot!(x_lin, runge.(x_lin))
scatter!(x, runge.(x), title = "Function approximation with piecewise constant polynomials", size=(800, 600))
function interp_chebyshev(x, xx)
vander_chebyshev(xx, length(x)) * inv(vander_chebyshev(x))
end
function interp_monomial(x, xx)
vander(xx, length(x)) * inv(vander(x))
end
function interp_error(ieval, x, xx, test)
"""Compute norm of interpolation error for function test
using method interp_and_eval from points x to points xx.
"""
A = ieval(x, xx)
y = test.(x)
yy = test.(xx)
norm(A * y - yy, Inf)
end
function plot_convergence(ievals, ptspaces; xscale=:log10, yscale=:log10, maxpts=40)
"""Plot convergence rates for an interpolation scheme applied
to a set of tests.
"""
xx = LinRange(-1, 1, 100)
ns = 2:maxpts
fig = plot(title="Convergence",
xlabel="Number of points",
ylabel="Interpolation error",
xscale=xscale,
yscale=yscale,
legend=:bottomleft,
size=(1000, 700))
for ieval in ievals
for ptspace in ptspaces
for test in [runge]
errors = [interp_error(ieval, ptspace(-1, 1, n), xx, test)
for n in ns]
plot!(ns, errors, marker=:circle, label="$ieval/$ptspace")
end
end
end
for k in [1, 2, 3]
plot!(ns, ns .^ (-1.0*k), color=:black, label="\$n^{-$k}\$")
end
fig
end
plot_convergence (generic function with 1 method)
So how accurate is it?#
plot_convergence([interp_monomial, interp_chebyshev, interp_nearest], [LinRange, CosRange])
Note that this convergence plot is specifically for the Runge function example given above, i.e., \(f(x) = \frac{1}{1+x^2}\). In general, convergence rates depend on smoothness of \(f\).
In general, a smoother f \(\Rightarrow\) faster convergence.