2025-10-13 Spline interpolation

2025-10-13 Spline interpolation#

Last time#

  • Runge phenomenon as ill conditioning

  • Stability and Chebyshev polynomials

  • Piecewise methods

  • Accuracy of piecewise constant (nearest neighbor) interpolation

Today#

  • Piecewise polynomial methods

  • Splines

  • Libraries

  • Generalization: boundary value problems

using LinearAlgebra
using Plots
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

runge(x) = 1 / (1 + 10*x^2)

CosRange(a, b, n) = (a + b)/2 .+ (b - a)/2 * cos.(LinRange(-pi, 0, n))

vcond(mat, points, nmax) = [cond(mat(points(-1, 1, n))) for n in 2:nmax]

plot_vcond(mat, points) = plot!([
        cond(mat(points(-1, 1, n)))
        for n in 2:30], label="$mat/$points", marker=:auto, yscale=:log10)

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] = x .* T[:, k-1]
        T[:, k] = 2 * x .* T[:,k-1] - T[:, k-2]
    end
    T
end

vander_chebyshev (generic function with 2 methods)

Piecewise 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

x = LinRange(-1, 1, 10)
s = LinRange(-1, 1, 100)
A = interp_nearest(x, s)
plot(s, A[:, 1:5])
scatter!(x, ones(length(x)))

I don’t see any visual artifacts#

x = LinRange(-1, 1, 30)
A = interp_nearest(x, s)
cond(A)

1.414213562373095

plot(s, A * runge.(x))
plot!(s, runge.(s))
scatter!(x, runge.(x))

function interp_chebyshev(x, xx)
    vander_chebyshev(xx, length(x)) / 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=(1200, 800))
    for ieval in ievals
        for ptspace in ptspaces
            for test in [runge]
                try
                    errors = [interp_error(ieval, ptspace(-1, 1, n), xx, test)
                             for n in ns]
                    plot!(ns, errors, marker=:circle, label="$ieval/$ptspace")
                catch
                    continue
                end
            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 maybe it’s not that accurate?#

plot_convergence([interp_monomial, interp_chebyshev, interp_nearest], [LinRange, CosRange])
plot!(xscale=:linear)

Observations#

  • Piecewise constant interpolation is well conditioned on any set of points

  • Piecewise constant interpolation converges very slowly (needs many points to increase accuracy)

  • Chebyshev/polynomial interpolation requires special input points, otherwise it is ill conditioned

  • Chebyshev/polynomial interpolation has “exponential” convergence for smooth enough functions

How could we improve the accuracy?#

Cubic splines#

  • Piecewise cubic function

  • Continuous values

  • Continuous derivatives

Splines#

If we are given an arbitrary distribution of points, interpolation with a single polynomial is not robust. Piecewise constant interpolation is not very accurate and gives a rough function. We could improve the accuracy by using a piecewise linear function, but the accuracy is still limited and the function still isn’t smooth (there is a “corner” where the slope changes at each data point). Splines are a way to guarantee an arbitrary amount of smoothness. The idea is that given sorted input points \(\{x_i\}_{i=0}^n\), we compute an interpolating polynomial \(s_i(x)\) on every interval \((x_i, x_{i+1})\).

Interpolation#

Given a function value \(y_i\) at each \(x_i\), we require

\[\begin{split}\begin{split} s_i(x_i) &= y_i \\ s_i(x_{i+1}) &= y_{i+1} \end{split} \end{split}\]
so that the polynomial interpolates our data. If the polynomial has order greater than 1, we are left with some extra degrees of freedom. To provide a unique solution, we’ll need to add conditions.

Smoothness#

The conditions above guarantee continuity, but not smoothness. We use our extra degree(s) of freedom to impose smoothness conditions of the form

\[\begin{split}\begin{split} s_i'(x_{i+1}) &= s_{i+1}'(x_{i+1}) \\ s_i''(x_{i+1}) &= s_{i+1}''(x_{i+1}) . \end{split}\end{split}\]
These conditions, which are applied at the interior nodes (\(x=1,\dotsc,n-1\)) couple the splines from adjacent intervals and causes the spline approximation to be globally coupled.

End-point conditions#

The conditions above are still not enough to guarantee a unique spline. Suppose we use quadratic polynomials for each \(s_i\). Then with \(n\) intervals, we have \(n\) degrees of freedom after imposing the interpolation condition. Meanwhile, there are only \(n-1\) internal nodes. If we impose continuity of the first derivative, we have \(n - (n-1) = 1\) undetermined degrees of freedom. We could fix this by imposing a boundary condition, such as that the slope at one end-point (e.g., \(s_0'(x_0)\)) was equal to a known value. This is not symmetric and is often an unnatural condition.

Suppose we use cubic polynomials. Now we have two degrees of freedom per interval after imposing the interpolation condition. If we impose continuity of the first and second derivatives, we have \(2n - 2(n-1) = 2\) remaining degrees of freedom. A common choice here is the “natural spline”, \(s_0''(x_0) = 0\) and \(s_n''(x_n) = 0\). Cubic splines are the most popular spline in this family.

Solving spline interpolation problems#

We need to choose a basis for the polynomials \(s_i(x)\). We could choose

\[ s_i(x) = a_i + b_i x + c_i x^2 + d_i x^3 \]
but this would be very ill-conditioned when the interval \((x_i,x_{i+1})\) is far from zero. A better-conditioned choice is
\[ s_i(x) = a_i + b_i(x - x_i) + c_i(x - x_i)^2 + d_i(x - x_i)^3 . \]
The interpolation property gives
\[\begin{split}\begin{split} a_i &= y_i \\ a_i + b_i(x_{i+1} - x_i) + c_i(x_{i+1} - x_i)^2 + d_i(x_{i+1} - x_i)^3 &= y_{i+1} \end{split}\end{split}\]
and continuity of the first and second derivatives gives
\[\begin{split}\begin{split} s_i'(x_{i+1}) = b_i + 2c_i(x_{i+1} - x_i) + 3d_i(x_{i+1}-x_i)^2 &= b_{i+1} = s_{i+1}'(x_{i+1}) \\ 2c_i + 6d_i(x_{i+1} - x_i) &= 2 c_{i+1} . \end{split}\end{split}\]
After trivially eliminating the \(a_i\), this is a block bidiagonal system (\(3\times 3\) blocks). We can reduce this to a scalar tridiagonal system. Define \(\delta_i = x_{i+1} - x_i\) and \(\Delta_i = y_{i+1} - y_i\). Then eliminate \(d_i\) using
\[ d_i = \frac{c_{i+1} - c_i}{3\delta_i} \]
and \(b_i\) using
\[\begin{split}\begin{split} b_i\delta i &= \Delta_i - c_i\delta_i^2 - \underbrace{\frac{c_{i+1} - c_i}{3\delta_i}}_{d_i} \delta_i^3 \\ b_i &= \frac{\Delta_i}{\delta_i} - \frac{\delta_i}{3}(c_{i+1} + 2c_i) . \end{split}\end{split}\]
Substituting into the equation for continuity of the first derivative gives
\[ \frac{\Delta_i}{\delta i} - \frac{\delta_i}{3}(c_{i+1} + 2c_i) + 2c_i\delta_i + (c_{i+1} - c_i)\delta_i = \frac{\Delta_{i+1}}{\delta_{i+1}} - \frac{\delta_{i+1}}{3}(c_{i+2} + 2c_{i+1}) \]
which reduces to
\[ \delta_i c_i + 2(\delta_i + \delta_{i+1}) c_{i+1} + \delta_{i+1} c_{i+2} = 3\left(\frac{\Delta_{i+1}}{\delta_{i+1}} - \frac{\Delta_i}{\delta_i}\right) . \]
To impose boundary conditions, we add a dummy interval on the right end (the actual value of \(x_{n+1}>x_n\) cancels out) so that the equation above is valid for \(i = 0,\dotsc,n-2\) and the right boundary condition \(s_{n-1}''(x_n) = s_n''(x_n)\) becomes \(c_n = 0\). The left boundary condition \(s_0''(x_0) = 0\) yields \(c_0 = 0\), so we must solve

(13)#\[\begin{equation} \begin{bmatrix} 1 & & & & & & \\ \delta_0 & 2(\delta_0+\delta_1) & \delta_1 & & & & \\ & \delta_1 & 2(\delta_1+\delta_2) & \delta_2 & & & \\ & & \ddots & \ddots & \ddots & & \\ & & & & \delta_{n-2} & 2(\delta_{n-2}+\delta_{n-1}) & \delta_{n-1} \\ & & & & & & 1 \\ \end{bmatrix} \begin{bmatrix} c_0 \\ c_1 \\ \vdots \\ \\ c_{n-1} \\ c_n \end{bmatrix} = \begin{bmatrix} 0 \\ 3\left(\frac{\Delta_1}{\delta_1} - \frac{\Delta_0}{\delta_0} \right) \\ \vdots \\ 3\left(\frac{\Delta_{n-1}}{\delta_{n-1}} - \frac{\Delta_{n-2}}{\delta_{n-2}}\right) \\ 0 \end{bmatrix} . \end{equation}\]

After solving this equation for \(c_i\), we will recover \(b_i\) and \(d_i\), and then can evaluate the spline at arbitrary points.

function spline_interp_and_eval(x, s)
    n = length(x) - 1
    function s_interp(y)
        delta = x[2:end] - x[1:end-1] # diff(x)
        Delta = diff(y)
        T = zeros(n+1, n+1)
        T[1,1] = 1
        for i in 2:n
            T[i, i-1:i+1] = [delta[i-1], 2*(delta[i-1] + delta[i]), delta[i]]
        end
        T[end,end] = 1
        rhs = zeros(n+1)
        rhs[2:end-1] = 3*(Delta[2:end] ./ delta[2:end] - Delta[1:end-1] ./ delta[1:end-1])
        c = T \ rhs
        S = zeros(n, 5)
        S[:, 1] = x[1:end-1]
        S[:, 3] = c[1:end-1]
        S[:, 5] = y[1:end-1]
        S[:, 2] = diff(c) ./ (3 * delta)
        S[:, 4] = Delta ./ delta - delta/3 .* (2*c[1:end-1] + c[2:end])
        S
    end
    
    function polyval(p, x)
        f = p[1]
        for c in p[2:end]
            f = f * x + c
        end
        f
    end

    function s_eval(S, s)
        f = zero(s)
        for (i, t) in enumerate(s)
            left = max(1, searchsortedfirst(S[:,1], t) - 1)
            f[i] = polyval(S[left, 2:end], t - S[left, 1])
        end
        f
    end
    
    A = zeros(length(s), length(x))
    aye = diagm(ones(length(x)))
    for i in 1:length(x)
        e = aye[:, i]
        S = s_interp(e)
        A[:, i] = s_eval(S, s)
    end
    A
end

# Cubic spline interpolation

spline_interp_and_eval (generic function with 1 method)

How well do splines work?#

x = LinRange(-1, 1, 14)
y = runge.(x)
s = LinRange(-1, 1, 100)
A = spline_interp_and_eval(x, s)
plot(s, [runge.(s) A * y])
scatter!(x, y)

@show cond(A)
plot(s, A[:,8])

cond(A) = 1.8495159246421367

Interpolations.jl#

using Interpolations

x = LinRange(-1, 1, 14)
y = runge.(x)
flin = LinearInterpolation(x, y)
fspline = CubicSplineInterpolation(x, y)
plot([runge, t -> fspline(t)], xlims=(-1, 1))
scatter!(x, y)

xx = LinRange(-1, 1, 100)
norm(runge.(xx) - fspline.(xx))

0.022048279334709518

function interp_spline(x, s)
    m, n = length(s), length(x)
    A = diagm(m, n, ones(n))
    for j in 1:n
        fspline = CubicSplineInterpolation(x, A[1:n,j])
        A[:,j] = fspline.(s)
    end
    A
end
s = LinRange(-1, 1, 200)
A = interp_spline(x, s)
@show cond(A)
plot(s, A[:, 1:5], legend=:none)

cond(A) = 1.962849489282676

Spline conditioning#

function my_spy(A)
    cmax = norm(vec(A), Inf)
    s = max(1, ceil(120 / size(A, 1)))
    spy(A, marker=(:square, s), c=:diverging_rainbow_bgymr_45_85_c67_n256, clims=(-cmax, cmax))
end

A = interp_spline(LinRange(-1, 1, 40), s)
cond(A)

1.7491451562213967

my_spy(A)
my_spy(abs.(A) .> 1e-6)

Accuracy#

plot_convergence([interp_monomial, interp_chebyshev, interp_nearest, interp_spline],
  [LinRange, CosRange], maxpts=60)
plot!(legend=:none)

Generalizations of interpolation#

To create a Vandermonde matrix, we choose a family of functions \(\phi_j(x)\) and a set of points \(x_i\), then create the matrix

\[ V_{ij} = \phi_j(x_i) .\]

Integrals?#

\[ B_{ij} = \int_{(x_{i-1} + x_i)/2}^{(x_i + x_{i+1})/2} \phi_j(s) ds\]

This leads to conservative reconstruction, which is an important part of finite volume methods, which are industry standard for shock dynamics.

Derivatives?#

What if we instead computed derivatives?

\[ A_{ij} = \phi_j'(x_i) \]
function diff_monomial(x)
    n = length(x)
    A = zeros(n, n)
    A[:,2] = one.(x) 
    for j in 3:n
        A[:,j] = A[:,j-1] .* x * (j - 1) / (j - 2)
    end
    A
end

diff_monomial(LinRange(-1, 1, 4))

4×4 Matrix{Float64}:
 0.0  1.0  -2.0       3.0
 0.0  1.0  -0.666667  0.333333
 0.0  1.0   0.666667  0.333333
 0.0  1.0   2.0       3.0

We need boundary conditions!#

First, a stable basis!#

Derivatives of Chebyshev polynomials also satisfy a recurrence.

function chebdiff(x, n=nothing)
    T = vander_chebyshev(x, n)
    m, n = size(T)
    dT = zero(T)
    dT[:,2:3] = [one.(x) 4*x]
    for j in 3:n-1
        dT[:,j+1] = j * (2 * T[:,j] + dT[:,j-1] / (j-2))
    end
    ddT = zero(T)
    ddT[:,3] .= 4
    for j in 3:n-1
        ddT[:,j+1] = j * (2 * dT[:,j] + ddT[:,j-1] / (j-2))
    end
    T, dT, ddT
end

chebdiff (generic function with 2 methods)

x = CosRange(-1, 1, 7)
T, dT, ddT = chebdiff(x)
c = T \ cos.(3x)
scatter(x, dT * c)
plot!(s -> -3sin(3s))

Solving a BVP with Chebyshev collocation#

A boundary value problem (BVP) asks to find a function \(u(x)\) satisfying an equation like

\[ -u_{xx}(x) = f(x) \]
subject to boundary conditions \(u(-1) = a\) and \(u'(1) = b\).

We’ll use the “method of manufactured solutions”: choose \(u(x) = \tanh(2x)\) and solve with the corresponding \(f(x)\). In practice, \(f(x)\) comes from the physics and you need to solve for \(u(x)\).

function poisson_cheb(n, rhsfunc, leftbc=(0, zero), rightbc=(0, zero))
    x = CosRange(-1, 1, n)
    T, dT, ddT = chebdiff(x)
    L = -ddT
    rhs = rhsfunc.(x)
    for (index, deriv, func) in
            [(1, leftbc...), (n, rightbc...)]
        L[index,:] = (T, dT)[deriv+1][index,:]
        rhs[index] = func(x[index])
    end
    x, L / T, rhs
end

poisson_cheb (generic function with 3 methods)

manufactured(x) = tanh(2x)
d_manufactured(x) = 2*cosh(2x)^-2
mdd_manufactured(x) = 8 * tanh(2x) / cosh(2x)^2
x, A, rhs = poisson_cheb(11, mdd_manufactured,
    (0, manufactured), (1, d_manufactured))
plot(x, A \ rhs, marker=:auto)
plot!(manufactured, legend=:bottomright)

“spectral” (exponential) convergence¶#

function poisson_error(n)
    x, A, rhs = poisson_cheb(n, mdd_manufactured, (0, manufactured), (1, d_manufactured))
    u = A \ rhs
    norm(u - manufactured.(x), Inf)
end

ns = 3:20
ps = [1 2 3]
plot(ns, abs.(poisson_error.(ns)), marker=:auto, yscale=:log10, xlabel="# points", ylabel="error")
plot!([n -> n^-p for p in ps], label=map(p -> "\$n^{-$p}\$", ps), size=(1000, 600))
Contents