25) Chebyshev Polynomials and Piecewise Interpolation

25) 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:#

  1. Chebyshev polynomials
    1.1 Chebyshev polynomials are well-conditioned

  2. Introduction 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))

ArgumentError: Package LaTeXStrings not found in current path.
- Run `import Pkg; Pkg.add("LaTeXStrings")` to install the LaTeXStrings package.

Stacktrace:
 [1] macro expansion
   @ ./loading.jl:1842 [inlined]
 [2] macro expansion
   @ ./lock.jl:267 [inlined]
 [3] __require(into::Module, mod::Symbol)
   @ Base ./loading.jl:1823
 [4] #invoke_in_world#3
   @ ./essentials.jl:926 [inlined]
 [5] invoke_in_world
   @ ./essentials.jl:923 [inlined]
 [6] require(into::Module, mod::Symbol)
   @ Base ./loading.jl:1816

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

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))

cond(ourvander(x_cos))

Definition of Chebyshev polynomials#

Define

\[ T_n(x) = \cos (n \arccos(x)) .\]
This turns out to be a polynomial, but it’s not obvious why. It is a trigonometric polynomial.

Recall

\[ \cos(a + b) = \cos a \cos b - \sin a \sin b .\]
Let \(y = \arccos x\) and check
\[\begin{split} \begin{split} T_{n+1}(x) &= \cos \big( (n+1) y \big) = \cos ny \cos y - \sin ny \sin y \\ T_{n-1}(x) &= \cos \big( (n-1) y \big) = \cos ny \cos y + \sin ny \sin y \end{split}\end{split}\]
Adding these together produces
\[ T_{n+1}(x) + T_{n-1}(x) = 2\cos ny \cos y = 2 x \cos ny = 2 x T_n(x) \]
which yields a convenient recurrence:
\[\begin{split}\begin{split} T_0(x) &= 1 \\ T_1(x) &= x \\ T_{n+1}(x) &= 2 x T_n(x) - T_{n-1}(x) \end{split}\end{split}\]

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

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))

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)

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

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.

Contents