2022-03-09 Higher dimensions
Contents
2022-03-09 Higher dimensions¶
Last time¶
Accuracy of piecewise constant (nearest neighbor) interpolation
Piecewise polynomial methods
Splines
Libraries
Today¶
Discussion
Compare accuracy and conditioning of splines
Generalization: boundary value problems
Cost for interpolation in higher dimensions
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
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
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
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]
vcond (generic function with 1 method)
Spline bases¶
using Interpolations
x = LinRange(-1, 1, 9)
y = runge.(x)
flin = LinearInterpolation(x, y)
fspline = CubicSplineInterpolation(x, y)
plot([runge, t -> fspline(t)], xlims=(-1, 1))
scatter!(x, y)
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, 100)
A = interp_spline(x, s)
plot(s, A, legend=:none)
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.5075717485062632
my_spy(A)
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=(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)
Accuracy¶
plot_convergence([interp_monomial, interp_chebyshev, interp_nearest, interp_spline], [LinRange, CosRange], maxpts=60)
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
Integrals?¶
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?
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, 10)
T, dT, ddT = chebdiff(x)
c = T \ cos.(3x)
scatter(x, dT * c)
plot!(x -> -3sin(3x))
Solving a BVP with Chebyshev collocation¶
A boundary value problem (BVP) asks to find a function \(u(x)\) satisfying an equation like
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(12, 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))
Curse of Dimensionality¶
Suppose we use a naive Vandermonde matrix to interpolate \(n\) data points in an \(n\)-dimensional space of functions, e.g., predicting \(z(x, y)\) from data \((x_i, y_i, z_i)\)
# A grid with 10 data points in each of d dimensions.
points(d) = 10. ^ d
flops(n) = n ^ 3
joules(flops) = flops / 20e9 # 20 GF/joule for best hardware today
scatter(1:10, d -> joules(flops(points(d))), xlims=(0, 10), yscale=:log10, legend=:none)
barrels_of_oil(flops) = joules(flops) / 6e9
scatter(1:10, d -> barrels_of_oil(flops(points(d))), xlims=(0, 10), yscale=:log10)
Fourier series and tensor product structure¶
For periodic data on the interval \([-\pi, \pi)\), we can use a basis \(\{ 1, \sin x, \cos x, \sin 2x, \cos 2x, \dotsc\}\), which is equivalent to \(\{ 1, e^{ix}, e^{i2x}, \dotsc \}\) with suitable complex coefficients. If we’re given equally spaced points on the interval, the Vandermonde matrix \(V\) (with suitable scaling) is unitary (like orthogonal for complex matrices) and can be applied in \(O(n \log n)\) (with small constants) using the Fast Fourier Transform. This also works for Chebyshev polynomials sampled on CosRange points.
points(d) = 10. ^ d
flops(n) = 5n * log2(n)
joules(flops) = flops / 20e9 # 20 GF/joule for best hardware today
scatter(1:10, d -> joules(flops(points(d))), xlims=(0, 10), yscale=:log10)
Partial differential equations¶
Boundary value problems in multiple dimensions.
Nonlinear (Newton), exploiting structure in linear algebra (matrix-free p-multigrid).
Lower-degree polynomials to fit noise-free data¶
We can fit \(m\) data points using an \(n < m\) dimensional space of functions. This involves solving a least squares problem for the coefficients \( \min_c \lVert V c - y \rVert \)
function chebyshev_regress_eval(x, xx, n)
V = vander_chebyshev(x, n)
@show cond(V)
vander_chebyshev(xx, n) / V
end
ndata, nbasis = 30, 20
x = LinRange(-1, 1, ndata)
xx = LinRange(-1, 1, 500)
C = chebyshev_regress_eval(x, xx, nbasis)
plot(xx, [runge.(xx), C * runge.(x)])
scatter!(x, runge)
cond(V) = 30.08350663794034
S = svdvals(C)
scatter(S, yscale=:log10)
