2025-10-15 Spline Interpolation¶
Accuracy of piecewise constant interpolation
Piecewise polynomial methods
Splines
Libraries
See the FNC
[1]:
using LinearAlgebra
using Plots
using Polynomials
default(lw=4, ms=5, legendfontsize=12, xtickfontsize=12, ytickfontsize=12)
# Here's our Vandermonde matrix again
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
# And our "bad" function
runge(x) = 1 / (1 + 10*x^2)
# And a utility for points distributed via cos
CosRange(a, b, n) = (a + b)/2 .+ (b - a)/2 * cos.(LinRange(-pi, 0, n))
# And a helper for looking at conditioning
vcond(mat, points, nmax) = [cond(mat(points(-1, 1, n))) for n in 2:nmax]
[1]:
vcond (generic function with 1 method)
Chebyshev interpolating polynomials¶
[2]:
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
[2]:
vander_chebyshev (generic function with 2 methods)
[3]:
x = LinRange(-1, 1, 50)
plot(x, vander_chebyshev(x, 5), title="Chebyshev polynomials", label=["\$T_0\$" "\$T_1\$" "\$T_2\$" "\$T_3\$" "\$T_4\$"])
[3]:
Well-conditioned¶
The Chebyshev polynomials are well-conditioned.
[4]:
a = 20
plot(range(2, a), [vcond(vander, LinRange, a)], yscale=:log10, legend=:topleft, xticks=range(2, a, step=2), xlabel="rows", ylabel="cond(A)", label="Vander, LinRange")
plot!(range(2, a), [vcond(vander, CosRange, a)], yscale=:log10, label="Vander, CosRange")
plot!(range(2, a), [vcond(vander_chebyshev, LinRange, a)], yscale=:log10, label="Chebyshev, LinRange")
plot!(range(2, a), [vcond(vander_chebyshev, CosRange, a)], yscale=:log10, label="Chebyshev, CosRange")
[4]:
Piecewise polynomials¶
Let’s take a different approach. We’ll re-examine the Lagrange polynomials as we vary the points.
[5]:
x = CosRange(-1, 1, 11) # What about LinRange?
s = LinRange(-1, 1, 100)
A = vander_chebyshev(s, length(x)) /
vander_chebyshev(x)
plot(s, A[:,1:5], size=(1000, 600), label=["\$l_0\$" "\$l_1\$" "\$l_2\$" "\$l_3\$" "\$l_4\$"])
scatter!(x[1:5], ones(5), color=:grey, label="values")
scatter!(x, zero(x), color=:black, label="roots")
[5]:
Let’s look at a specific # Fit the function |x| y = abs.(x) plot(s, [abs.(s) A * y], label=[“true” “fit”])example that seems like it should be easy to fit but shows strange behavior.
[6]:
# Fit the function |x|
y = abs.(x)
plot(s, [abs.(s) A * y], label=["true" "fit"])
[6]:
The effects of the higher order functions seem to be producing different behavior between the points we are fitting than we want. The absolute value function looks like two linear functions glued together. Let’s try that.
New interpolation matrix¶
[7]:
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))
[7]:
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
[8]:
A = interp_nearest(x, s)
plot(s, A[:, 1:5], label=["\$c_0\$" "\$c_1\$" "\$c_2\$" "\$c_3\$" "\$c_4\$"])
scatter!(x, ones(length(x)), label="values")
[8]:
Visually the accuracy seems reasonable.
[9]:
x = CosRange(-1, 1, 41)
A = interp_nearest(x, s)
plot(s, A * runge.(x), label="fit")
plot!(s, runge.(s), label="true")
scatter!(x, runge.(x), label="data")
[9]:
This looks pretty good, but how accurate is it really?
[10]:
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
[10]:
plot_convergence (generic function with 1 method)
plot_convergence([interp_monomial, interp_chebyshev, interp_nearest], [LinRange, CosRange])
Piecewise interpolation seems to be not especially accurate, but it doesn’t have the runaway Runge phenomenon.
Observations¶
Piecewise constant interpolation is well-conditioned on any set of points
Piecewise constant interpolation converges very slowly (needs many points)
Chebyshev/polynomial interpolation requires special input points otherwise it is ill-conditioned
Chebyshev/polynomial interpolation has “exponential” convergence for ‘smooth enough’ functions
Cubic 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). High(er)-order piecewise polynomial interpolation can be very useful, such as for finite element methods, but we will take a different direction here.
Splines are a way to guarantee an arbitrary amount of smoothness. The idea is that given sorted input points \(\left\lbrace x_i \right\rbrace_{i = 0}^n\), we compute an interpolating polynomial \(s_i \left( x \right)\) on every interval \(\left( x_i, x_{i + 1} \right)\).
This approach gives us a
piecewise cubic function
continuous values
continuous derivatives
Interpolation¶
Given a function value \(y_i\) at each \(x_i\), we require
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
These conditions, which are applied at the interior nodes (\(x = 1, \dots, 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 - \left( n - 1 \right) = 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' \left( x_0 \right)\)) 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 \left( n - 1 \right) = 2\) remaining degrees of freedom. A common choice here is the “natural spline”, \(s_0'' \left( x_0 \right) = 0\) and \(s_n'' \left( x_n \right) = 0\). Cubic splines are the most popular spline in this family.
Solving problems with spline interpolation¶
We need to choose a basis for the polynomials \(s_i \left( x \right)\). We could use
but this choice would be very ill-conditioned when the interval \(\left( x_i, x_{i + 1} \right)\) is far from zero.
A better-conditioned choice is
$
The interpolation conditions give us
and continuity of the first and second derivatives gives
We can trivially eliminate the \(a_i\), as \(a_i = y_i\). This leaves 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
and \(b_i\) using
Substituting into the equation for continuity of the first derivative gives
which reduces to
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, \dots, n - 2\) and the right boundary condition \(s_{n - 1}'' \left( x_n \right) = s_n'' \left( x_n \right)\) becomes \(c_n = 0\). The left boundary condition \(s_0'' \left( x_n \right) = 0\) yields \(c_0 = 0\), so we must solve
After solving this equation for \(c_i\), we will recover \(b_i\) and \(d_i\), and then can evaluate the spline at arbitrary points.
[11]:
# Cubic spline interpolation
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
[11]:
spline_interp_and_eval (generic function with 1 method)
Let’s try it¶
[12]:
# Problem to fit
x = LinRange(-1, 1, 14)
y = runge.(x)
s = LinRange(-1, 1, 100)
# And the fit
A = spline_interp_and_eval(x, s)
plot(s, [runge.(s) A * y], label=["true" "fit"])
scatter!(x, y, label="data")
[12]:
Our new approach is well-conditioned and low error.
[13]:
@show cond(A)
plot(s, A[:,8], label="\$s_7\$")
cond(A) = 1.8495159246421362
[13]:
Interpolations.jl¶
See the documentation here.
We can do this all with the Interpolations package instead of hand-coding our \(A\).
[14]:
using Pkg
pkg"add Interpolations"
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), label=["true" "fit"])
scatter!(x, y, label="data")
Resolving package versions...
No Changes to `~/.julia/environments/v1.11/Project.toml`
No Changes to `~/.julia/environments/v1.11/Manifest.toml`
[14]:
And the condition number matches up.
[15]:
# Make a matrix so we can compute condition number
s = LinRange(-1, 1, 100)
A = zeros(100, length(x))
for i in 1:length(x)
y = zero.(x)
y[i] = 1
f = CubicSplineInterpolation(x, y)
A[:, i] = f.(s)
end
@show cond(A)
plot(s, A[:, 8], label="\$s_7\$")
cond(A) = 1.849515924642137
[15]: